{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Linear regression baseline\n",
    "\n",
    "In this notebook, we will create the linear regression baselines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import seaborn as sns\n",
    "import pickle\n",
    "from src.score import *\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "sns.set_style('darkgrid')\n",
    "sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def to_pickle(obj, fn):\n",
    "    with open(fn, 'wb') as f:\n",
    "        pickle.dump(obj, f)\n",
    "def read_pickle(fn):\n",
    "    with open(fn, 'rb') as f:\n",
    "        return pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Load and prepare data for training\n",
    "\n",
    "First up, we need to load and prepare the data so that we can feed it into our linear regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "DATADIR = '/data/weather-benchmark/5.625deg/'\n",
    "PREDDIR = '/data/weather-benchmark/predictions/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Load the entire dataset\n",
    "z500 = xr.open_mfdataset(f'{DATADIR}geopotential_500/*.nc', combine='by_coords').z\n",
    "t850 = xr.open_mfdataset(f'{DATADIR}temperature_850/*.nc', combine='by_coords').t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Load the validation subset of the data: 2017 and 2018\n",
    "z500_valid = load_test_data(f'{DATADIR}geopotential_500', 'z')\n",
    "t850_valid = load_test_data(f'{DATADIR}temperature_850', 't')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Split into train and test data\n",
    "z500_train = z500.sel(time=slice('1979', '2016'))\n",
    "z500_test = z500.sel(time=slice('2017', '2018'))\n",
    "t850_train = t850.sel(time=slice('1979', '2016'))\n",
    "t850_test = t850.sel(time=slice('2017', '2018'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Compute normalization statistics\n",
    "z500_mean = z500_train.mean().values\n",
    "z500_std = z500_train.std('time').mean().values\n",
    "t850_mean = t850_train.mean().values\n",
    "t850_std = t850_train.std('time').mean().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array(54111.594, dtype=float32),\n",
       " array(1115.3833, dtype=float32),\n",
       " array(274.54495, dtype=float32),\n",
       " array(5.681409, dtype=float32))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z500_mean, z500_std, t850_mean, t850_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# Normalize datasets\n",
    "data_z500_train = (z500_train - z500_mean) / z500_std\n",
    "data_z500_test = (z500_test - z500_mean) / z500_std\n",
    "data_t850_train = (t850_train - t850_mean) / t850_std\n",
    "data_t850_test = (t850_test - t850_mean) / t850_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(32, 64)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, nlat, nlon = data_z500_train.shape; nlat, nlon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f77c48e1c50>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEcCAYAAADtODJSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2df3wU9Z3/XzOTbH5vQsBgAtRUvmhjEeXEYj3QyhWhGsTW60mJevbwrEXQs3oa9QqItRq1Wm0B+wPto23a4zyvKCko9MR6Vy2VlquU1B+lUDBZiEDI7x/7Y75/hGzzec/szG52N7vZvJ6Pxzwe+exn5jPv+cxk3zvzec3ro5mmaYIQQggZJnqqAyCEEDK6YSIhhBASF0wkhBBC4oKJhBBCSFwwkRBCCIkLJhJCCCFxwURCoqK5uRkzZ85EMBhMdSiEkDSDiYTYMm/ePLzxxhvhckVFBfbs2QPDMFIYlTM//vGP8bnPfQ7Tp09HbW2tpf7555/H/PnzMXPmTCxbtgxHjx4N1910002YOXNmeJk+fToWLVoUrv/d736Hv//7v8fMmTOxaNEi7N692zGWDz74ANdffz3OO+88LFy4UOnLlpYW3HLLLZgzZw7OPvtsfPDBB67HtmXLFlx22WU4//zzsXz5cpw8eTJct3XrVixZsgTnnXcerr/+ete2CEk0TCQkYygrK8Py5ctxzTXXWOp+85vf4IknnsD69euxa9cuTJ48GXfeeWe4/vvf/z727NkTXmbOnIkFCxYAAE6ePIkvf/nLWLZsGXbv3o2bbroJX/7yl9HW1hYxljvvvBPnnHMOdu3ahTvuuAO33XYbTpw4AQDQdR1z587Ft771raiO6/3338eqVavw6KOP4le/+hXy8vLwwAMPhOtLSkpwww034J//+Z+jao+QRMNEQiz867/+K5qbm3HLLbdg5syZ+N73vocPPvgAZ599NgKBAADg+uuvx5NPPoklS5Zg5syZuOWWW9Da2oo777wTf/M3f4NrrrlG+aW9f/9+fPGLX8QnPvEJLFiwAFu3bk143Jdffjk+/elPo6SkxFK3c+dOLFy4ENOmTYPH48Hy5cvx1ltv4dChQ5Z1P/jgA+zevRuLFy8GAOzZswcTJkzAZz7zGRiGgcWLF6O0tBTbt2+3jePAgQPYt28fVq5cidzcXCxYsABnnXUWXnnlFQDAhAkTUFNTg3PPPTeq49qyZQvmzZuHCy+8EAUFBbj99tuxY8cOdHZ2AgAuvvhiXHHFFZg4cWJU7RGSaJhIiIXHHnsMFRUVeOaZZ7Bnz56Iv3S3bt2KRx99FK+//joOHTqEJUuW4JprrsFvfvMbTJ06FevWrQMAdHd345/+6Z9QXV2NN954A0888QQeeOABvP/++7btrlmzBrNmzbJdhj5uigXTNGHnBvTee+9ZPtu8eTNmzZqFKVOmRNzWNM2I8f/pT3/ClClTUFhYGP7sYx/7GP70pz8NK/b3338fZ599drj8kY98BNnZ2Th48OCw2iMk0TCRkGHzuc99Dh/5yEdQVFSESy65BFOmTMHFF1+MrKwsLFy4EI2NjQCA1157DZMmTcI111yDrKwsfPzjH8eCBQvCv9Ala9aswe7du22XLVu2DCvWSy+9FNu2bcM777yD3t5erFu3Dpqmobe317Luiy++iM9+9rPh8syZM9HS0oKGhgb4/X787Gc/w6FDh2y3BYCuri4UFRUpnxUVFaGrq2tYsXd3d1vaKywsHHZ7hCSarFQHQEYvEyZMCP+dk5OjlHNzc9Hd3Q0AaGpqwttvv41Zs2aF64PBIK666qoRi/WTn/wkbrvtNtx2223o6OjAjTfeiIKCApx++unKert378axY8fC4yMAMG7cOKxfvx51dXVYu3Yt5syZg4svvjj8KOnKK69Ec3MzAOB73/seCgoKwo+dBuns7ERBQYFrnLt37w7fAVZUVODnP/858vPzh90eISMBEwlJOuXl5bjwwgvx3HPPRbX+qlWrIt55DH65DoeamhrU1NQAGBjH2LBhA6ZNm6ass3nzZsyfP9/yJf2JT3wCL7zwAgAgEAhg/vz5+OIXvwgAlngOHDiAw4cPo7OzM/x465133kF1dbVrjLNmzcKePXuUz6ZNm4Z33nknXD58+DD8fj8qKyujOGpCkg8fbRFbJkyYgMOHDyekrU996lM4ePAgNm/eDL/fD7/fj7fffhv79++3XX/t2rWKgmro4pREAoEA+vr6EAqFEAwG0dfXFxYH9PX14b333oNpmmhubsaqVatwww03oLi4OLx9b28vXn75ZeWx1iCNjY3w+/3o7OxEXV0dJk6ciLlz59rG8dGPfhRVVVVYt24d+vr6sGPHDrz77rvKXU5fXx/6+/sBAP39/ejr64t4XIsWLcLOnTuxe/dudHd346mnnsL8+fPDSWrosYZCIfT19cHv90dsj5BEw0RCbLn55puxYcMGzJo1Cxs3boyrrcLCQmzcuBFbt27F3LlzMWfOHDz++OPhL9JEsWHDBsyYMQPf/e538dJLL2HGjBnYsGEDgIEv7jvvvBMzZ87E5z//eZx//vm4/fbble1/8YtfoKioCBdddJGl7e9///u46KKLcOmll+LDDz8MCwki8cQTT+APf/gDLrzwQjz++ON4+umnUVpaGq6fMWMGZs6cCQD4zGc+gxkzZkRsa9q0aXjggQdw11134eKLL0ZXVxdWr14drn/xxRcxY8aM8NjSjBkz8NWvftW9wwhJEBontiKEEBIPvCMhhBASF0wkhBBC4oKJhBBCSFwwkRBCCIkLJhJCCCFxwURCCCEJpuvESfeVMogxI/+96dldaGkfeOnLYsAXUtcNhUSXRNFFbqtommsTrhvouvqZJn4GaGIbwxAr6NY2DV1dx8gS+zCc28zOsv4W8Yh1PKLNbEsbhqgX+7TpC9mGR8Rh6HKfmmO93Wd6jCdNbh+U15ENIXHhyG1kOSSuVQDoD4bEOqZjvSz7g9Y4/YGgWEe0EXCOU65v1xemPDb5TxTr/6nNOkHZYSHZplouK8rBd//pE5Z2h8MTc/8erR8ciVg/bvLp+Mr//GdC9pVqxoxFSkt7H3xtAyZ78uKxXNCyPppEYvMPPhT5pe+GTAqAXSJxLssvfVlvt45MJLql3vkL3O6zHJdtPCKRWJKCTV+47UN+qbslGrvP7NZxYjiJxC1xuJUBoD/g/KXdF3BOJHL7gc+CjuvINuU+3WICYv+/c1t/YB21HBTH6va/n0jam46g7VBTxHoj1h+XacyYSSSEEDKS6Jpzsojxd0paw0RCCCFJwNA027vpofWZAhMJIYQkgWxdg8fhtiM7g25JmEgIISQJGC6PtjhGMgqpnFKMwtI8AEBHb0Cp6xbl3m7Vgru7Q7X47um0utb2dKiz1QV61ImIgv09ajmgtmFd3zr7ntzGjSxPnlI2PLmWdYwcdZ3s3EKlnJWnlj356jwdufkeS5u5BdmO6xQUqeUyb45SLsxRL8vxhdZ9FOer+ygS2+RnGzGVBz4TA/iGFAE4q8ksIjm4f1OEIAfT1fqAGAzu7FevVQBoE9dva496/bb1qfUfiuu5s9faZofYxk04YBE3SJGGza/vLBdxg5t4wS4m2V+WQX8xgN/Tr4oKJhRYr7XhYsDl0VYU18doIa3eI9m5cyeuvvpqLF68GIsWLcL27dsBDEwUdO2112LBggW49tprOVc1ISTt0TDwBRtpyZw0kkZ3JKZp4u6770Z9fT3OOussvPPOO/jCF76AT3/601i9ejWWLl2KxYsX48UXX8SqVavwwx/+MNUhE0JIRMbSYHta3ZHouo6Ojg4AQEdHB8rKytDa2orGxsbwNKXV1dVobGzEiRMnUhkqIYQ4MjhG4rRkCmlzR6JpGr75zW9i+fLlyM/PR1dXF77zne/A5/Nh4sSJME49rzYMA2VlZfD5fMqMc4QQkk6MJdVW2tyRBAIBfOc738H69euxc+dObNiwAXfccQe6u7tTHRohhMSMfurRVqQlVguedCZt7kj++Mc/oqWlBRdccAEA4IILLkBeXh5ycnJw9OhRBINBGIaBYDCIlpYWlJeXpzhiQgiJDOW/KeD000/HkSNH8Oc//xlnnnkm9u/fj2PHjuGMM85AVVUVGhoasHjxYjQ0NKCqqirmx1qzzxyPzlNSP2kO1+NXJYAnhLy3Rcgl/3JMlfoCQNsJVRZ78kO13HnsmFIOtH+olKXc1y/kwABghtQ4Qy5yYH9Xm2M9AOhZqtxRyoE9+V6lLOXAXUJiDAA5xacp5bwidZ3eblXu2yekpt1CkinPD2Aj6yyUpoHOvkx2vwZzhGQ1N0ctF3rUcl6WlLyqZSlvBaySVk36SWnOktdASJU9A0CPMFDs8qt909ar9l9Th3qt+cT1DQBdQmbsJsWV9bkuRpyA9ZxI70i/nUOlQwx2ccgBbXneZViFHqssfLiMpcH2tEkkp512GtasWYPbb789bFj48MMPo6SkBGvWrEFtbS3Wr18Pr9eLurq6FEdLCCHO6HDx2hqxSJJP2iQSALjqqqtw1VVXWT6fOnUqnn/++RRERAghw8Ojw3Gw3ZNBmSStEgkhhGQKdP8lhBASF8kYI+nr68PXv/51vPnmm8jJycH555+PBx98MJ4wEwITCSGEJIFk3JE89thjyMnJwSuvvAJN03BMiHhSBRMJIYQkgUTfkXR1dWHz5s345S9/GRYkTZgwIa4YE8WYSSTZhhaWIBpi3lu3N0ylzLCz129Zp6tDleLKaW2l+29/R6tSlnJfPdvqQqrrqjRRbiPlwZpYX5btPjOEHFjWe/KLlXJ+qSr1BQCPcOLNEpLKoJDu9oi+O9qv1tu507Z1q9sUd6hxnzE+X92gUJUcZ0mrXgC5Yjpeq/zUeXrfHE1dXwvayLNDzteaJudkNlS5r8306pCzHeeKn8F6ntr/edlq30wqsrpCS5dh6SDsJv+N5mU7Oa+7lHm39an7aBOuxtK5F7DKwmVchblZjuXsUOJGwKN9j8Tn8yEYVI/F6/XC61Wl94cPH0ZJSQm+/e1vY9euXSgoKMDtt9+OWbNmJSzm4TJmEgkhhIwkWZqObD1yYso69aOhpqYGTU3q3O4rVqzAypUrlc8CgQAOHz6Mc845B/fccw9+//vf45ZbbsGOHTtQWKi+3zXSMJEQQkgS0AwNmsMtyWBdfX297R2JpKKiAllZWWED2/POOw/jxo3DgQMHcO655yYw8thhIiGEkCSgGxp0h0QyWBet3VNpaSlmz56NX/3qV5gzZw4OHDiA48eP44wzzkhIvPHAREIIIUlA0zVoNmNxQ+tj5YEHHsB9992Huro6ZGVl4dFHH7W9exlpmEgIISQZuDzaGo5r45QpU/CjH/0ojqCSAxMJIYQkAV13ebSVQa+2j5lEEgqZYbdRv2mjoRyCR+gpi/NVCWZpryolBYDu8epgWUBIGf19Feo+ilT34mikulLeK919/b2qHFi6A4f8VjmqbFO6+5ZM/n9KecIk9TZ6YpnqcgwAJfmqFDcoZLQdQs570sZ9Vt3eer5kG6WF1nMyFOn43BewOsue6LHKuociJa1S4TouVzreWp16JQEhKe4X+t7jXWpMIThfu3Zt9Ao3YHkcsm8Aq1tvWYFVjj4Ui9uyaNNv59Qr+i9bPAYalyekz8Wqi3S7jQz/QyELly7ERUKK3iskyHlSSx0HRrYBIzuym7BT3WhjzCQSQggZSQbGSBxUW7wjIYQQ4sSA/NdhsD2DZrZiIiGEkCSguch/mUiSRCRnywMHDqC2thYnT55ESUkJ6urqUFlZmepwCSEkIpqmOT6+0jhDYnKI5Gy5evVqLF26FIsXL8aLL76IVatW4Yc//GGKoyWEkMgYHgOGw9S9TnWjjbRJJJGcLY8fP47GxkY899xzAIDq6mo8+OCDOHHiREzztgfNvxreHetUVUJt3ar6IyAUJnLu7ZJ8qxpnfGGJUj5jvKpm+vAjqtmh74Rq4tgj5ok3bVQukpA5LqZtCkus5nxVU9S4P1ZepJRPEwqsQmHIaGd4KRU6ncLwr6ldnTP81/uPK+WjH3Yr5Q7RVwCQLeLo6lL77//eU+218wvV4zirwvoS17SJqmLNoioSxyqVSh19qgJIzvFu95n8LpGmjHJeeIupI4A+oTzqFHEcFX0j+7/5pLV/pdnhRK967UwtVY0fy4VqLicr9l/bsr/l1Syvq76A9XruFmrJTmHsaFHmBZ3nhY8HDrangEjOlrm5uZg4cSIMY+A/zjAMlJWVwefzxZRICCFkJBkYIxkbg+1pM2vwUGfL//qv/8Jdd92FlStXoru7231jQghJMwZNG52WTCFt7kgiOVvm5ubi6NGjCAaDMAwDwWAQLS0tURudEUJIKtA1zfHt9WjmbBktpM0dyVBnSwBhZ8vKykpUVVWhoaEBANDQ0ICqqio+1iKEpDWaobsumULa3JEAkZ0t16xZg9raWqxfvx5erxd1dXWpDpUQQhwxsnUYNoKLofWZQlolkkjOllOnTsXzzz+fgogIIWSYuLzZPhz333QlrRJJMunsC6D9lAzVTe5bJOZxzhMazXKvVXYo5++WUtFjYp9Sguw2tzQA5AuTt2IhgZUmd6WinGfzCyhXXOjysa308wtGYcYXEFNpd4j5v2Vfybm337IxVJT0iP6T8ul+ITk++WGXYxkA/iykz+cJafTlHytTygUe577r7Lceh1ynQJgKegy1L7o0tb7bb23ziDj2vUdV887DJ5wFK3bX2iRhkFhZopZPK1CvrUJxbeWLsp0w3e17tFdooeWx6zYSY49oNFdca9IAVJp3ZjlMjRsrdP8lhBASF5ru8mY7EwkhhBAndEN3fI/EqW60wURCCCFJgG+2E0IIiQst24CeHfkrVuPEVoQQQpzQdZdHWwkc2E81mXMkhBCSRiTrhcRvf/vbOPvss/Hee+8lOOLhM2buSM4r94alfnaS1aFI6aK0Msi2+SVRmKN+VuAy97OMQcp/c22erUrZrJQ6amL+dS2kSmDtRJimod5eC+NYizRaxh2QdrUAejX1s6Ap5NNFqlOsbNPXprrTdnZb5+YO+NVL15Mn5NO6KtHOEuf09PGqey1gdWyeWqa6AU8Qrs+lYo72wihswY2Aemx6+4fqCqZ6AgqKVMlxm40yuldIWP0htVwqnI8nC2lv1WnqcQPA6ULe69XVc6D1tqvlgHqtmbrav8Fcq9uylPP2iPJwfHmlDNlSllJ+cS3KuerjwS1ZDCeR7Nu3D//3f/+HioqKeEJLOLwjIYSQJDAwsZUeeYnRa6u/vx9r167F6tWr025SrDFzR0IIISOJ7smG7rHOXTS0HgB8Ph+CQfVpgtfrhder3sU99dRTuOqqqzBlypTEBxsnTCSEEJIENE1zHFAfvKuoqalBU1OTUrdixQqsXLkyXN6zZw/27t2Lu+66KznBxgkTCSGEJIFox0jq6+tt70iG8tZbb+HPf/4z/u7v/g4AcOTIESxbtgwPP/ww5syZk+DIY4eJhBBCkkC0iSSauZVuvvlm3HzzzeHyvHnz8Mwzz+Css86KP9AEwERCCCFJQNMGBtWd6jOFMZNIZpblwjzlXNsbUhUPPQFVfioFEdLJN8vG2iAXzlJbU1w0mnDR1fo7RVl1dAUA9LrIe01nwaRpeKyfZavyxzyxTihblYoGxC+sPt0q/w2a0lFV7a8i4Vo8rVSVn447Vx2gbDlzvGUffcJiWDqpSom2lE5LV1jAKg2dWKD2TXmhGncB1HOk9VkdhSWaX8h/pYw2KGS24hwXFloHWotz1P465zRVtiyPozRPPc6SXKtsOavNp5SNjqNKOdSlxq3nqecwJGTLdvJf6cTbL8o54lorsEjsrf+Hek+rUtYC4v9I/B+acr4Q3QBgjXU4aC428vFMtfvqq68Oe9tkMGYSCSGEjCR6dhbgifwV62SfMtpIy3sr+ebmgQMHcO2112LBggW49tprcfDgwdQGSAghLji+Q6I7P/YabaTdkdi9ubl69WosXboUr7zyCpYuXYpVq1alMEJCCHFn8NFW5CW9XiqMh7RKJHZvbh4/fhyNjY2orq4GAFRXV6OxsREnTpxIZaiEEOKIprt4bfGOJDnYvbnp8/kwceJEGKc8oQzDQFlZGXw+X6RmCCEk5fDRVgoYfHNz6dKlqQ6FEELiRjcM1yVTSBvZQKQ3N++9914cPXoUwWAQhmEgGAyipaUlqpd4hqLt/e+wPLMwv0ipKypQ5X6akMSauuimLKt/jpT3hnKLlbIe6FP3EbSR97pgkTJKt18p/5VlG926PDZTHJsu6nWPKvPMyh9nadPIkY6rzr9X5KPiyhKrTFniYuAclnpHws70Trop5/W3KWW9XZWWQvRNsECVKbeZ6nUEAMVC+mx0tCjlUIf6yFbr71HK42wk3kUTKpVyb1DtPym/9vSpx5XVdMjSZqD5z0q5r+24Gpf4Na2PU+W+hrhO9KDVwTlLV79Ipev2uCz1WI2TH6gxtKmSZAAItor+7FEl2aZwRtY8qks0cguAi662tDsctGwDuoNqK5MmtkqbO5Kbb74Z//u//4tXX30Vr776Kk4//XRs3LgRV1xxBaqqqtDQ0AAAaGhoQFVVFUpLS1McMSGERGYsPdpKmzsSJ9asWYPa2lqsX78eXq8XdXV1qQ6JEEIc0XQdcHohkYkk+Qx9c3Pq1Kl4/vnnUxgNIYTEhmZozokkg+S/aZtICCFkNKPpOuDktcU7EkIIIU5oRjaQFVk4ohmRJ70abYyZRBI8ehhm94BaRctRlRpaboFzWayvS6UHAE0ov7K6TzrGY1WPqOaIoTxV9QUAWqBDbcOvKsEsOiWp6rJRm2niCtD6nNVkerdQLrUfsayTJY5FqpsgjAi1gLNRoSmUSwAAoWAz/UIVJI7dFGWErOonU6rihKoI4rrQiyeoTeaqasBiG1WO3qUqpuSxmb3djmVdHgeAbNFfRo4ah+YX+2h+Xyn3HT1saTMoVFpmUCioclRFmkX9JK5/w+Y6KRT9pfeofaO1qwqsQPMBNUah0AKAYJdqfhrySzNVFel3pRWUwKq1GyYudySOdaOMMZNICCFkJNEMA3B4V0TjeySEEEIc0XTrXa2szxCYSAghJBnohnMicaobZTCREEJIEqBqixBCSHy4qLZA1RYhhBBH3GxQeEcyCtEjD3yZ/eo82lI6KsuhgNWADt2qNDfmi0TEZifVDQrJqkWuKtDEryG7NiVSJmvtGxFD0CpH1bKd9xPqEjJmsQ+zT8h9o5DqSjm1lKtaYrR541j+00tJqy4k3oaIwZAGmFIGDavEVZoKyrnQTXGtyfUBQO9wlpoHhJQ3JMqBTlUyC1hls7K/LGVx/QePq3Jfw0a2LP8bpeRYxhkSx+nvtPZFqN85biNX/Z+wypjdDUOjxuE7J1yfIYydREIIISNJggfbW1tbcffdd+PQoUPweDw444wzsHbt2rQwsM2clEgIIWmEphvQDIclxkSiaRpuuukmvPLKK9iyZQumTJmCxx9/PEnRxwYTCSGEJANd++vb7bZLbKaNJSUlmD17drh8/vnno7m5OdFRDws+2iKEkCSgGR4gy8GixRgYj/H5fAiKsUav1wuv12u3FQAgFArhpz/9KebNm5eQWOOFiYQQQpLB4B2JUz2AmpoaNDU1KVUrVqzAypUrI2764IMPIj8/H9ddd11CQo2XtEkkTgNJBw4cQG1tLU6ePImSkhLU1dWhsrIy1SETQkhEtCgH2+vr623vSCJRV1eHv/zlL3jmmWegp4nySzPdJrceIU6ePIl33303/Aywrq4ObW1t+PrXv44bbrgB11xzDRYvXowXX3wRL7zwAn74wx/G1H7Ptu/C7D4lrZRyUin7lGZqUkYbxcmTclQLNnJIV+Q28iJ1q49icE8em5Sfmn4hObaRQksJsZSSuklzXfvOpo1Y5b7RyH/lfNt6br5azhNuwGLecl24BQNWqbNF7tsr5hiX/RvFdSOdkKVUOprz4XYOpGuukSf6pqhEKVvcgWFzbclj7xN9Jfuu39nZF7CR/4o4pVRdy/ci9/Jlru1GQ2DfLwE75+pBPHnI+vilMbX55JNP4ne/+x2++93vIi/PKi9PFemRzhB5IOn48eNobGxEdXU1AKC6uhqNjY04ceJEqkIlhBB3NJfBdi22wfb3338fzzzzDFpaWrBkyRIsXrwYt956a5KCj420ebQ1lKEDST6fDxMnToRx6i7BMAyUlZXB5/OlhX6aEELs0LI9ABzu7LJje/lx2rRpePfdd+MLKkmkzR3JUNJtIIkQQmJm8M32iEtafv0Oi7S7I5EDSeXl5Th69CiCwSAMw0AwGERLSwvKy8tTHSohhERE01xmSMyg+UjS6kiefPJJ/OEPf8C6devgOeV5M378eFRVVaGhoQEA0NDQgKqqKj7WIoSkN453Iy6KrlFG2tyRDA4kVVZWYsmSJQCAyZMnY926dVizZg1qa2uxfv16eL1e1NXVpThaQghxQdOc7zpiHGxPZ9JG/ptserash9k94LxqcbgVUkeLC+wwZLQxMxw5cCIQx+LmKCyRUlI7huPEG2+bbvuIyv3X5bxbXI6jcVt2Oc+ucms7pFQ3xmvJzs061v6VcmDpPB3VeICLu3U0snDLOXSaDwSw9JWWX4y8Rctd9xMNoUNvA07/T1ke6B+ZkZB9pZq0uSMhhJBMwtSzAN0h+emZ8/WbOUdCCCHpRgY9vnKCiYQQQpKBrgMmZ0gkhBAyTExNdxlsZyIhhBDiBBNJ5mGaQYtaaxA3E0a5ne1Tz3iVXG4GjMlCHpvLcVj6YhiKq+FsI7GohGKMe1hItZ8UOwkHV8vc8zZYDEJd6s2gzXHI6zfGRya6Xd9ZlI3O/Wfpf7l+Avo/GrNUC24qOVlvJvD/TjcAJ1Es3yMhhBDiBB9tEUIIiROXRJJexiJxwURCCCHJQNOc5b8ZJA1mIiGEkGQwhixSmEgIISQJDIyROAy2c4yEEEKII0aW810HVVujD02L3rbZInGNxrTRbb70RDASF16M88LbSTItJpguct/hyDpdDRVzxBzhYv5vO1mo1czT5ToYjsRVmmS6rC7lv3ZyYSkJdpMUWyTEdv0vzQ9DcUQs1WAAAB+NSURBVMrTo7l245SiDwe5D01L4P8YH20RQgiJDz3CS2dD6jOEqI9k27Zttp+//PLLCQuGEEIyBVPTXZdMIeojuf/++20/X7VqVcKCceLAgQO49tprsWDBAlx77bU4ePDgiOyXEEKGxeALiU5LjKTr96DrkRw+fBiHDx+GaZrhvweXN954IzwlbrJZvXo1li5dildeeQVLly4dsQRGCCHDQsNf3yWxXWJvMl2/B10Tyfz583H55Zejp6cH8+fPV5a7774bK1euTHqQx48fR2NjI6qrqwEA1dXVaGxsxIkTJ5K+b0IIGRZ6FkyHJdaJrZL1PXjDDTfg+eeft3x+8803R92G65G88847AIDrrrsOP/7xj2MIL3H4fD5MnDgRxik1imEYKCsrg8/nQ2lpaUpiIoQQRzQdjrq8U6otn8+HoFDSeb1eeL1e5bNkfQ/u2bMHx48fxx//+Efcf//94fZ3794ddRtRp8RUJZGRwE1GOCLuv9G0l4p53V0kmXZ9J+W8FjnwCEzoY0Yh95XE7HwcjSzcDZd5yN3kwbbbxCgHtiXOeeCt7aXg2sUwzlkCr01T05zP36lEUlNTg6amJqVqxYoVI/K0BwCys7OxadMm3HnnnbjxxhvxrW99CyUlJTCdnIsFUSeSQCCAn/zkJ3jrrbfQ2tqq7KS+vj62yGOkvLwcR48eRTAYhGEYCAaDaGlpQXl5eVL3Swghw8U0nV3kB6mvr7e9I5Ek83uwsLAQzzzzDJ544glcc801WLduHbQY3nOJOv0+/PDD2LRpE2bNmoV9+/bh8ssvx/Hjx3HRRRcNK/BYGD9+PKqqqtDQ0AAAaGhoQFVVFR9rEULSlpBpui7AQIKYPHmystglkmR9Dw7eFGiahjvvvBNf+cpX8MUvfhH9/f1Rt6GZUd6/zJ07F5s2bUJFRQVmzZqF3bt3Y//+/Vi9evWIPPbav38/amtr0d7eDq/Xi7q6Opx55plRb9+zZT3M7jbbuljfkLV9BJKKt87TgGj6LhmPtmJ+rJSAvkvKoy03EtFX0UyOJUn0o60kMJw3210dEfK9yLvylnjCCtPe3ev4eEjTNHjzcyPW2xHv96AdW7duxRVXXKF81tjYiFdffRUrVqyIqo2oH2319vaGb6Fyc3PR09ODqVOnorGxMYaQh8/UqVNtlQWEEJKWmKZzIhlGk8n4HpRJBADOOeccnHPOOVG3EXUimTp1Kvbu3YsZM2Zg+vTp+Na3voXCwkJMnDgx6p0RQshYIWQOLGOBqO+b77vvPmRlDeSd2tpaNDY24rXXXsODDz6YtOAIIWS0YkaxZAqOdyRvvvlmxM++9KUvAQD8fn8Swko8oaAfZmB4sVrlrDaSV7dGRsINeDjPseOMK5rbcy1LuB8kY7xCrpCV7dxAVGM7McrCh+Fwa8FyTsVYRTKk0y4SZDuGdexxYnVnjj1uC/J8BBP3fTaW7kgcE0kkf62haJqG//7v/05YQIQQkgmYLmMkmZRjHBPJq6++OlJxEEJIRhECEHTKFpkzHQnnIyGEkGTg9mjLaRbe0QYTCSGEJAHXR1tMJIQQQpwInVoikUFPtphICCEkGbh5bfGOZBQS8gdh9geiWlczhNw3qP6u0LOt3eZq15AAyaurxDIRcuBYiUb2KSSrWraQ5gp5sJRbh7o6LE3GKjeV+zSDNjLaGOXhUtas5RWoK9j0v9nTpZblOi7nLCp7HilXj8YSRSLjkFYibjFEc+25bOMm95X/l3bIdeT/tkTzJ+5/xhzip2WHzkRCCCHEiaDprNriYDshhBBH+GiLEEJIXIRgIuTw2mEC3stPG5hICCEkGbhNbMU7EkIIIU6E4PxCIu9IEswDDzyAN998Ex6PB/n5+bj//vtx7rnnAgB6enpw7733Yt++fTAMA/fccw8uu+yymPcR6vMj1Dsw45dUblhUWi4TMYUQnforFtzUJIA1Lj3GszdiZpMWlZYnwor2SHNNLcdm8h+3uALOs7vZzWNu+Z/vdzYmlEowzWIUaWMcKc5BsF2dbM1NiaR7rCdd98RpmGhnBCkVU6I/ZV9Y6qMx6oxTVWinntQ86rUS6u1WysFeNU6LIjMr+lkB3QiaJoJUbY0cl1xyCe677z5kZ2dj586duOOOO/CLX/wCALBx40YUFBRgx44dOHjwIGpqarB9+3YUFBS4tEoIIaljLA22J8GTOnYuu+wyZJ/6dXf++efjyJEjCJ369b1t2zYsWbIEAFBZWYnp06fj9ddfT1mshBASDabpPG97JiWStLgjGUp9fT0+9alPQT91u93c3IxJkyaF68vLy3HkyJFUhUcIIVERDA0sEeszyCNlRBLJZz/7WTQ3N9vWvfHGGzBOPa/++c9/ji1btqC+vn4kwiKEkKQRcnmz3alutDEiieRnP/uZ6zo7duzAk08+iR/84AeYMGFC+POKigo0NTWhtLQUAODz+TB79uykxUoIIYkg5DLYnuhE4iRaSjZpMUayc+dOPPzww9i4cSMmT56s1C1cuBCbNm0CABw8eBB79+7F3LlzUxEmIYREjT9kwh90WBI8D+8ll1yCLVu24KWXXsKXvvQl3HHHHQlt34m0GCO59957kZ2djdtuuy382Q9+8AOMGzcOy5YtQ21tLebPnw9d17F27VoUFhbGtT+L/DcJ82DHbBYnJbNSPglYZZpyzmoppxRyVM1OjuqG27zYdsaEQgpq9vfGtEuLGaKN/FdKbaXs07LPKOb71j1S5i1wMXXUC7zqPmzMEkNd7WIddS8hvyotlxJX3WYueosJppQpu81fb4OUYMPFbNJyvcpr1VZi7DwfvTYMlb0ujDOl/Ff2d7BfPU7Tnzhpf7SPtnw+H4LiWvF6vfB6vXabRWToaxFDRUt6Er7fJGmRSH79619HrMvPz8fTTz89gtEQQkj8RPtoq6amBk1NTUrdihUrsHLlymHvW4qWkk1aJBJCCMk03KbaHayrr6+3vSORpLNoiYmEEEKSQDBkIuiQSQbrysvLo2ovHtFSsmEiIYSQJGDCeYzETLBr46Bo6bnnnrOIlpINEwkhhCQBfwiOyqysBLs2OomWkg0TCSGEJIGQy6OtUILlv06ipWQzZhKJZ1wxzJwBTwIpJ5XySIs0VMpZ7ebAdpPJuiDn+9bzi6wrScfV3i7rOkPblLJPG0mxnQuusg+/s5TXNga/kFSK/pMSV0tM0uXYce3BldTjkFJcKYW2d6OVclS1zVDbcaVsiuMMinrbfbhJiIXc18h1mRcegJarfial0LbuyW6IOM38HrXsJumORoLs5hAsyqHuDjWGHuu1J8+Bxd1X9G9WgeirBJrB8s12QgghceE2Z7tT3WiDiYQQQpKA6XJHYvKOhBBCiBODFikR6w0mEkIIIQ5wjIQQQkhcjLT7byphIiGEkCQQCpmOEt9Ey39TyZhJJFmTpwL9AxJGizxSSGDNPiFxlVLIKNxsLRJhF2mjlEJaHF1t0LJLlLIuZaBCKqrl5Fvj1NVLQAup0lyzT7inyr6x6QuLRFj0hZszr6WsWyXKUvJqFIm+KB6v7jNLSGA1m+npxC9EvVWdiTMg5agdJ0W5VW3ORhJukaN6nN19NSFj1nOt59AiHXe5DmxlyRKbPndqQ0qh5fWr5Rdbm8gRcZmib/pVyXGo44RablfLABASkmBDSM8t/2fiuwC5iZP/UrVFCCEkLkw4P77KoDzCREIIIcmgPxhCv8Ok7U51o420mCFxkF27dqGqqgo//vGPw5/19PTgX/7lXzB//nwsXLgQO3fuTGGEhBASHYMWKZEWjpEkgc7OTjz++OO45JJLlM83btyIgoIC7NixAwcPHkRNTQ22b9+OggRaGRBCSKKJ1kY+E0ibO5JHHnkEy5YtszhVbtu2DUuWLAEAVFZWYvr06Xj99ddTESIhhERNyHS5I8kg+W9aJJJf/vKXaG9vx8KFCy11zc3NmDRpUrhcXl6OI0eOWNYjhJB0wimJuN2tjDZG5NGW0xSRL7/8Mr7xjW/gueeeS2oM5hkzwpLFkEd9LGYawoVUyBC1oJA2+lVJ7MBnfWrZUq9uo3Wq0kWzq10ph2xktdJB1RhXppSDhacpZX/JJKXc0WeVfUoJYqFH/W3hCagSTL1XjVMLqMcNWI8dmsvvFSE5ln0pzwcAi5tvyFOolAO5RaJeSGL96nEBgN7XqZZFvSHPiZCSmt1q39i5REupuUWKLqW70rnXTgotP7Nz2nXYJwzr+lIW7ioHNtT15f9UMNfqZt3jUSXB8os1X1iI6ONU91+9Vy0DgNEnHIHltSOdvqUs3EjcV6I/YKI/EHlA3R9gIokJpykid+/ejQ8//BCf//znAQCtra3YuXMnTp48iRUrVqCiogJNTU0oLS0FAPh8PsyePXskwiaEkGETNF3GSDLo0VbKB9tnzZqFN998M1yura3F9OnTcd111wEAFi5ciE2bNuHcc8/FwYMHsXfvXnzjG99IVbiEEBIVY2mwPeWJxI1ly5ahtrYW8+fPh67rWLt2LQoLC903JISQFDLSMySmkrRLJI888ohSzs/Px9NPP52iaAghZHgEQiYCDsnCqW60kXaJhBBCMoGQyxhJJsl/x0wiORIqCP8CaGpRjdyOdauKKrcTXJxrVblMKVYft+VlqZqf3FxVy1RcrCqqjI6jSjnLRpEilTQBr6raaoFq8LfvL6qK6EinVWElB/yKc9R9TPKqqpaJBaoZYm6+VZHVL6Rg7f2qeqmjT1VplYj+LMhV28zLshosasJ0sduvqmN6etWyR8wTXySVSwC8eUJRJRRpxoQKpawXqEokqdLScvIs+9Cyc9RtskQ5W5pLuiv0tYCqJnNTvVkUWTYqL6lmCuWpCiuLCk6afQrVVmufVb109KQa54keVR1piHN8epFqWOnNsyrB8orUbeR1IpHXqq5pGB9h3VhJlUXKrl27cOONN+L+++8PjzUnmzGTSAghZCQJhZwH1G3MoeMmkkNIsmEiIYSQJBCt/Nfn8yEo7ma9Xi+8Xq/dZo4MOoS89tprMW8bD0wkhBCSBKJVbdXU1KCpqUmpW7FiBVauXBnT/oY6hDCREEJIBhAMhRB0eH41WFdfX297RyJJB4eQSDCREEJIEoj2hcTy8vKo2ovHISTZMJEQQkgS6A+G0OfgtZVI1ZabQ0iyGTOJZHdzB7pOyVAPtapy3+OdqhxY/orIEVLe0kJh8gigqV2VS47LUyWV44TE9bQCtVxWNEUpFxTamDYKGWenrsohT3So8sk+caH62q1tdvaqss0ppWqb8h+hrVc99vxsq5lfa68aR5PYb7auSjKzDSH3FW1KSTIA+MU56uxXj0Puo7xIldkG863n0BT7KSxRz4lRpMqtJUFD3UePzZdIj5Apd4l1ZL30fSzwWOXAJV7ZX2rZkEab/cJ01EbubkqZcrYqZe4Nqf3bH1T7rldIvttsDEOPdavXSZu4btqETPz946oho7xuAPf/O7lJr+j/3CwdFxdZZdvDIRWqrVQxZhIJIYSMJKn02pIOIcmGiYQQQpIA3X8JIYTEBU0bCSGExEUgEELAYbDdqW60wURCCCFJIGSajncdNG0khBDijGnCdEoWTCSJ50c/+hHq6+uRnZ0NwzCwefNmAEBPTw/uvfde7Nu3D4Zh4J577sFll10Wc/t/PtaFk6fcRbuFNLHHr5Y7hQzRd1KVr3qyrLLDoly1K4tt5KVD+Vi56lz66akTlPJp+ar8EgAMYWTaJSSVfuFkWuhRZaAzK1QHV8D9V5Gc090QslopMQas0tvTC9VjkdLdgGhDSnnbxfmwwyIFFedIzleva9Y2pYS1VWxjmupxdfWrcbd0tanb91j30SGuvWPCkVnKseW1KqXoAHCaV+3fM4WEe6ooj8tVy3buylniHEoZsnRbbhVx68J119dhdZ4+3q3K7mXfyPPeJuTCduMPhZb/Q+f562V/F+dm4eLKxPj/miHAdLgjkdPJj2bSIpFs374dL7/8Mv7zP/8ThYWF+PDDD8N1GzduREFBAXbs2IGDBw+ipqYG27dvR0FBgUOLhBCSWkIhl0dbGTTY7j7ZwQjw7LPPYsWKFeEpdE877bRw3bZt27BkyRIAQGVlJaZPn47XX389JXESQki0DNyROC+ZQlokkv379+P3v/89lixZgs997nP4j//4j3Bdc3MzJk366yRQ5eXlOHLkSCrCJISQqAmGQggGHZYMerV9RB5tOblWvvHGGwgGg/D5fPjJT36C1tZWfOELX8BHP/pRXHjhhSMRHiGEJJ6Q6ThGggx6tDUiicTJtRIAKioqUF1dDV3XMX78eFx88cV4++23ceGFF6KiogJNTU0oLS0FMDAJzOzZs0cibEIIGTamSyJxTDKjjLR4tFVdXY3/+Z//AQB0d3fjt7/9LT72sY8BABYuXIhNmzYBAA4ePIi9e/di7ty5KYuVEEKiIWSarkumkBaqrRtvvBFf/epXceWVVwIAFi9ejL/9278FACxbtgy1tbWYP38+dF3H2rVrw4Pyw8VOQjmUfvHGaZmQV9rJDqUk2CPkqJNLVUdRKYltEzJEO4Si0iJD7/Y7y4GzpX4YQH62egnkuvSNbLOj3xq3lPfKsnQUlnH3BtSyXX9LGbLU8MlzLGXKdtLczn6rQ+1QOoQb7REh3X3H16GUT3RZJa+SlnZ1nZYm1am3v0fdZ2GJ6jINAOPL1KP/YKIqLf/LhB51fSFN99q4K5dKF11Rll+C1utC7V95DQCAVzjzynMmr1d5zqV0F7D+737Y7nwOZJty+3gwTZc7EiaSxJKbm4vHHnvMti4/Px9PP/30CEdECCHxwfdICCGExEXolDrLqT5TYCIhhJAkYJrOdx0Z9GSLiYQQQpLBWHqznYmEEEKSwFiS/zKREEJIEmAiyUCKcrMwaN4q5aSyPL5QlUdKiaAs2yHlwEVCYpkl5MHS8dZOYy7lkHKszi8sF2Qb3X5rm1KKW+hR47STDA/FTi4spaBS3tvWp0pDpUSzQ8g67eTaecLZuN+jHkev6Jxs3VlaGs068rj6xT6s14k1bmmLYZExe1V5b26Bur4nRz1uwHqtSTfrE52qy6683v02g75SriudeV0uCwt214l0iS4S57RYyIPLhWReXlcDn8UnLS+wkUIPl1DIZbCdFimEEEKcSIX8N9J0HMmGiYQQQpLASA+2O03HkWyYSAghJAmYLjMkDtb5fD4ExcxhXq8XXq83pv09++yzuP32222n40g2TCSEEJIEorVIqampQVNTk1K3YsUKrFy5Mqb9DU7H8dRTT6G/vx9LlizBP/zDP8Qe+DBgIiGEkCQQ7aOt+vp62zsSSTpPxzFmEklhTlbY9dDNdVOqSXSpxpHuiTbbSFWWrHdTfvXamMf5QyIu0YSMSyqw5LTmdsi5tiWy7+zitJvHfSi5IpBSoZKTiqweGzNF+Zk08Au4KPPskOekSMz/LecDl2opacwpzT7t4iwvVrcxpjhfJ3Kfdp/J/pMGohZzRBt1WUCcw25RL1Vb8rrJFvvUbfpfnhOLEaSLuixoc0ql0q5IU/siP1stS6VjoceqihsuZsCPUKDfoX5AlVZeXh5Ve/FMx5Fs0sJGnhBCMg3TDMIMOSyms9t0rDhNx5FsmEgIISQJmMEQzGDQYUms/vfGG2+Ez+fDlVdeic9//vNYtGhReDqOZDNmHm0RQshIMnjn4VSfSJym40g2TCSEEJIERjqRpJK0SCQHDhzAqlWr0N7ejv7+flxxxRVh6VtPTw/uvfde7Nu3D4Zh4J577sFll12W4ogJIcQZ0ww5J5IMmtkqLRLJY489hgULFuC6665DV1cXqqurcemll2LGjBnYuHEjCgoKsGPHDhw8eBA1NTXYvn07Cgrk5KqEEJI+mIF+F9VW5LrRRlokEk3T0NExMN91b28vNE1DaWkpAGDbtm145JFHAACVlZWYPn06Xn/9dXzmM5+JaR85WbpFEjqIVbooylIObCORlTJNt22s+1B1D3ZSXTdprp0seSjBKGbSkRJMN2NIaWQ4sI3zfnKynCWWUq4q5ayAdW5tWZbyYDejTsB6DuU6djLkoUgZbrGY5zzSZ07Ic25nlugmT3drQ24f1TYxysSlmaIdUu4r25DvZNjJf+U2bnHK/7ssGyn0cAmFQgg53JFkkmljWqi27rvvPmzduhVz587FvHnzsGzZMkyePBkA0NzcjEmTJoXXLS8vx5EjR1IVKiGERMVIy39TyYjckbi9kblp0yYsXrwYN910E1paWnD99ddj+vTpOO+880YiPEIISThmyGWMJIPuSEYkkbi9kfmjH/0Iv/jFLwAAZWVluOiii/DWW2/hvPPOQ0VFBZqamsKPunw+H2bPnp30mAkhJB4G3xdxqs8U0uLR1uTJk8NvZHZ2duK3v/0tpk2bBgBYuHAhNm3aBAA4ePAg9u7di7lz56YsVkIIiYZQcMAiJeIS9Ls3MkpIi8H2hx9+GF/72tfw7LPPIhAI4IorrsCll14KAFi2bBlqa2sxf/586LqOtWvXhm2SCSEkbXF5jwR8jySxTJ8+Hf/+7/9uW5efn4+nn356hCMihJD4GBgjiTwOwjGSUUhxbrZFVhoJd5lt/G1I9anFMTcBP1aicbx1c0K2k1jGsj1glULLsnSjdZN92uEWpySauN1wO8fREOvc53ZS3VgZTtyyf92k5FLKGw2yLyxy9mEce6zHKt2B44EvJBJCCIkLWqQQQgiJi1Ao6PJCIhMJIYQQB0IBP0L+yDYooQBVW4QQQpxweSERHGwffeTaTFEaCbcxvWgGSDW4DNgnYNDUjWgG203EO9hu85mLZ5LbPkw52B7FwPiYGWxPwD61BLSRDGLti2iItb9i+Z5wo+I0r2MiqTjNOi/7aEUz5X8tIYQQEgNp8WY7IYSQ0QsTCSGEkLhgIiGEEBIXTCSEEELigomEEEJIXDCREEIIiQsmEkIIIXHBREIIISQumEgIIYTExZiwSDlw4ABqa2tx8uRJlJSUoK6uDpWVlSMex7x58+DxeJCTkwMAuOuuuzB37twRja+urg6vvPIKmpqasGXLFpx11lkAnPso2fFFiilSf41ETK2trbj77rtx6NAheDwenHHGGVi7di1KS0tT1ldOMaWyrwBg+fLl+OCDD6DrOvLz8/HVr34VVVVVKb2uIsWU6r7KSMwxwPXXX29u3rzZNE3T3Lx5s3n99denJI7LLrvMfPfddy2fj2R8b731ltnc3GyJxSmGZMcXKaZI/TUSMbW2tpq//vWvw+VHHnnEvPfee133ncy4nGJKZV+Zpmm2t7eH/96xY4d59dVXu+472XFFiinVfZWJZHwiOXbsmHnBBReYgUDANE3TDAQC5gUXXGAeP358xGOxu4BTFd/QWJxiGMn4ok0kqeizl19+2fzHf/zHtOmroTGZZnr11c9+9jPzs5/9bFr11WBMpplefZUpZPyjLZ/Ph4kTJ8IwBqbQNAwDZWVl8Pl8KC0tHfF47rrrLpimiQsuuABf+cpX0iI+pxhM00xpfLK/vF7viPdZKBTCT3/6U8ybNy9t+mpoTIOkuq/uv/9+/OpXv4Jpmvj+97+fFn0lYxok1X2VaXCwfQSpr6/HSy+9hBdeeAGmaWLt2rWpDimtSZf+evDBB5Gfn4/rrrsuJfu3Q8aUDn310EMP4bXXXsMdd9yBRx99dMT3b4ddTOnQV5lGxieS8vJyHD16FMHgwLwAwWAQLS0tKC8vT0ksAODxeLB06VL87ne/S4v4nGJIZXx2/eUWb6Kpq6vDX/7yF3zzm9+Erutp0VcyJiA9+mqQq6++Grt27cLpp5+e8r6SMbW2tqZVX2UKGZ9Ixo8fj6qqKjQ0NAAAGhoaUFVVNeK3qt3d3ejo6AAwMHHT1q1bUVVVlRbxOcWQqvgi9ZdbvInkySefxB/+8AesW7cOHo/Hdd8jEZddTKnuq66uLvh8vnD51VdfRXFxcUr7KlJMOTk5Kb+uMpExMbHV/v37UVtbi/b2dni9XtTV1eHMM88c0RgOHz6MlStXIhgMIhQKYerUqfi3f/s3lJWVjWh8X/va17B9+3YcO3YM48aNQ0lJCX7+8587xpDs+OxieuaZZyL210jE9P7776O6uhqVlZXIzc0FAEyePBnr1q1LWV9Fiqm2tjalfXXs2DEsX74cPT090HUdxcXFuOeee/Dxj388ZX0VKSav15vSvspUxkQiIYQQkjwy/tEWIYSQ5MJEQgghJC6YSAghhMQFEwkhhJC4YCIhhBASF0wkJGOZN28e3njjjVSHQUjGw0RCCCEkLphICCGExAUTCcl4+vv78dBDD2HOnDmYM2cOHnroIfT39wMAdu3ahUsuuQTPPvssPvnJT2LOnDl44YUXUhwxIaMLJhKS8WzYsAG///3v8eKLL+Kll17C3r17sX79+nD9sWPH0NHRgddffx0PPfQQ1q5di7a2thRGTMjogomEZDxbtmzBrbfeivHjx6O0tBS33norXnrppXB9VlYWbr31VmRnZ+PSSy9Ffn4+Dhw4kMKICRldMJGQjKelpQUVFRXhckVFBVpaWsLlkpISZGX9dY63vLw8dHd3j2iMhIxmmEhIxlNWVobm5uZw2efzhd1eCSHxw0RCMp4rr7wSGzZswIkTJ3DixAmsW7cOixYtSnVYhGQMGT9nOyHLly9HV1cXrrrqKgDAwoULsXz58hRHRUjmwPlICCGExAUfbRFCCIkLJhJCCCFxwURCCCEkLphICCGExAUTCSGEkLhgIiGEEBIXTCSEEELigomEEEJIXDCREEIIiYv/D2WJUMf5CMc3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_z500_train.isel(time=0).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f77a00ac748>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEcCAYAAADtODJSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2de3wU5bnHfzOzl2QTNiGBhIRbCiJGQeCApa0CQkWiBuKlLQiirXgrBq09HA3aIkJVUrVaFLAq4lGix2OtgEEFPFClVRGtVRRRwEQgWYgkIQm5bPYy54+QNe8zk5lNNpvdbJ7v57OfT959Z2eeeWd2n8y8v/k9kqqqKhiGYRimk8iRDoBhGIbp2XAiYRiGYUKCEwnDMAwTEpxIGIZhmJDgRMIwDMOEBCcShmEYJiQ4kUSYadOm4b333gvrNgoKCvDoo4+GdRuR5KOPPsKMGTMiHQbD9Fo4kTCdQlVVPProo5g0aRLGjx+P+fPn48CBA4H++fPnY/To0Rg3bhzGjRun+aF///33kZOTgzFjxmD+/PkoKysLetsjR47Et99+G2hPmDABW7duDX2nwshjjz2GmTNn4uyzz8bjjz8u9KmqirVr1+LCCy/Ef/zHf+COO+7AqVOnAv2XXXZZYBzHjRuHs88+G7fcckugf8eOHcjNzcW4ceMwZ84cHDx40DCWL7/8EldeeSXGjBmDK6+8El9++WWg7+uvv8aCBQswceJEjBw5Mqh9e+6553D++edj/PjxWLJkCZqbmwN9GzZswJVXXolRo0ahoKAgqPUxPQ9OJEynePPNN/Hqq6/ixRdfxIcffoixY8fizjvvFJZZunQpPvnkE3zyySfCD31VVRXy8/Nx++2348MPP8SoUaNwxx13dPcudCtDhw7F4sWLMWXKFE3fxo0bsWnTJrz00kvYtWsXmpqasGLFikD/li1bAuP4r3/9CxkZGcjJyQEAlJaWYvHixVi2bBn27NmDqVOn4te//jW8Xq9uHM3NzVi4cCFmzZqFPXv24PLLL8fChQsDP/4WiwU5OTm4//77g9qvXbt24amnnsJzzz2HHTt24OjRo1i1alWgPy0tDQsXLsRVV10V9FgxPQ9OJFGE3+/HU089hYsuuggTJ07E7bffjpMnTwIAFixYgA0bNgjLz5o1C9u2bQMAHDp0CL/61a/wwx/+EDNmzMAbb7wR1liPHj2K8ePHY/DgwVAUBbNmzTL9T7iV7du3Y8SIEbjkkktgt9uxaNEi7N+/H4cOHTL97Lx58wAAeXl5GDduHN544w3s3r0bkydPDiwzbdo0PPPMM5g5cybGjh2Lu+++GydOnMANN9yAcePG4Ze//CVqamoCy//73//GnDlzMGHCBMyaNQu7d+/u4GiYc8UVV2DKlClISEjQ9O3cuRM/+9nPkJGRgYSEBNx4441444030NjYqFl2z549qKqqwsUXXwwA+Mc//oEJEyZgwoQJsFgsuPHGG3H8+HHs2bNHN44PP/wQXq8X1113HWw2G6699lqoqooPPvgAADBs2DD8/Oc/x4gRI4Lar40bN+JnP/sZRowYgaSkJCxcuBCvvfZaoP/iiy/GRRddhOTk5KDWx/RMOJFEEc8//zzefvttbNiwAbt27UJSUhKWL18OAJg5cyaKi4sDyx48eBDl5eW48MIL0dDQgOuvvx65ubl477338Kc//Qn33XefcKupPT766KPAD5He66OPPtL93GWXXYbDhw+jpKQEHo8Hr732GiZNmiQs88gjj2DixImYM2eO8ON84MAB4baJw+HAkCFDgkpERUVFAIBNmzbhk08+waWXXqq73LZt27B+/Xps3boVO3fuxI033ojf/va32L17N/x+P1544QUAwPHjx3HzzTfj17/+NT788EPcdddduO2221BVVaW73ptvvrndsbr55ptN49dDVVW0dSpSVRXNzc3C7btWXnvtNcyYMQMOh6Pdz6qqiq+//lp3WwcPHsTIkSMhSVLgvZEjRwb9TwDlwIEDOOuss4R1nThxAtXV1Z1aH9MzsUQ6AOZ7Xn75ZSxduhQDBgwAAOTn52Pq1Knwer246KKLsGzZMpSVlWHgwIF4/fXXMX36dNhsNrz99tsYOHBg4PbBOeecgxkzZmDr1q2m/1kaJQsj+vfvj/HjxyMnJweKomDAgAH47//+70D/4sWLMXz4cNhsNmzZsgW33HILNm3ahCFDhqChoQEpKSnC+hITE1FfX9/hONrjmmuuQb9+/QC07GNKSgrOPvtsAMD06dPx/vvvA2hJSJMnTw7ccjr//PMxatQovPPOO7jiiis06/3LX/7SZTG2MnnyZDzzzDO45JJLkJSUhKeffhoANFckjY2N2Lp1K9auXRt47yc/+QkeeeQR7N69G+PGjcPTTz8Nj8eDpqYm3W3V19ejT58+wnuhjH1DQwMSExMD7dZ119fXo2/fvp1aJ9Pz4EQSRZSXl+PWW2+FLH9/oSjLMiorK5Geno4pU6Zgy5YtuOmmm7Bly5bAffSysjJ89tlnmDBhQuBzPp8Ps2bNClusq1evxueff4533nkH/fr1w+bNm3Hddddhy5YtiI+Px5gxYwLLXnHFFSguLsY777yD+fPnw+FwCJPJQMsPj95tn87SmkQAwG63C+24uDg0NDQAaBnzt956Czt37gz0e71eTJw4sctiMeOqq66Cy+XCtddeC6/Xi+uvvx47d+4M/EPRyrZt25CcnIwf/vCHgfeGDx+OlStXYsWKFfjuu+8wc+ZMnHHGGUhPTwcAjBs3LrDsli1bkJCQ0Omx37x5M+69914AwPjx4/HMM89ojmXr3115LJnohxNJFDFgwAA88MADGD9+vG5/bm4unnjiCZx33nloamoK/NhlZGTgvPPOw/r16zu8zY8++gg33nhju/1PP/20kKBa2b9/Py655JLAj92VV16JBx54AAcPHsTo0aM1y0uSFLgFM2LECOE+ekNDAw4fPowzzjijw/GHSkZGBvLy8vCHP/whqOVvuOEGfPzxx7p9rT+uHUWWZdx222247bbbALTMe6SnpweSQSsbN25EXl6ecFsKAHJycgKT77W1tXj11VcDx+CTTz4Rlj3jjDPw7LPPQlXVwHq++uorzJ071zTOWbNmaf45GTFiBL766qvALcb9+/ejX79+fDXSy+A5kiji6quvxmOPPRaQwlZVVeHtt98O9E+ZMgXl5eVYtWoVLr300sCVy4UXXojS0lJs3LgRHo8HHo8Hn332WVCT1xMmTAgogvReekkEAEaPHo233noLJ06cgN/vx8aNG+H1ejF06FDU1tZi165dcLvd8Hq92Lx5Mz766CNccMEFAFpuLR04cABbt26F2+3G6tWrMXLkSAwfPhwA8Le//Q3Tpk1rN+Z+/frhyJEjwQ2qCbNmzcLOnTuxa9cu+Hw+uN1u7N69G8eOHdNd/plnnml3rIySiMfjgdvthqqq8Hq9cLvd8Pl8AICTJ0/i8OHDUFUVBw8exMqVKzVXpseOHcPu3bt1b7d9/vnn8Pl8qKqqwtKlSzFt2rTAWFJ++MMfQlEUPP/882hubg4IOH70ox8BaJljcbvd8Hg8AAC32y3IeSl5eXn461//ioMHD6KmpgZr164VYmzdV7/fHxjf9hRlTM+Fr0iiiFYFzfXXX4+Kigqkpqbi0ksvxUUXXQQAsNlsmD59Ol599VVBLpuYmIh169Zh5cqVWLlyJVRVxciRI7FkyZKwxXrjjTeisrISl19+ORoaGjB06FCsWrUKTqcTVVVVeOyxx/DNN99AURQMGzYMq1evxrBhwwAAKSkpePzxx7F8+XL813/9F8aMGYM//elPgXW7XC7hlgwlPz8fBQUFaGpqwvLly5Gamtrp/cjIyMCaNWvw0EMP4T//8z8hyzLOPfdcLFu2rNPr1OP3v/+9cBX25JNP4sEHH8SVV16J6upq3HLLLTh27BhSUlJw7bXXYvbs2cLnN23ahLFjx2LIkCGadd9///3Yv38/rFYrcnJyDJ/XsNlsWL16NX73u9/hkUcewfDhw7F69WrYbDYALbdJf/rTnwaWP/fcczFw4EDs2LFDd32TJ0/GDTfcgGuvvRZNTU2YMWNG4MoKANauXYsnnngi0N68eTPy8/OxaNEikxFjehISF7Zioo3rr78e99xzT7v/VTMME11wImEYhmFCgudIGIZhmJDgRMIwDMOEBCcShmEYJiQ4kTAMwzAhwYmEYRimi6mvOhnpELqVXqPauuqBbThW3eJdpKo+oU/1+cW232fY9vs82g3QZfxknWSbXYEkKUK77QNsnUZWDLsl0i9J2m0q9jiySqvYtoiPL9G4ZVl8cltnE9rPWEzWITZ13ggDOl8t+pbPS889cQE//YDeOv2at8Rt0PObrMOv84CganL+0nOPItFzUee8oscoHPjJeNLvKWVA33j8teCiLtn2nyb9DNVH9R9sBYC+gwbgt7v+2iXbijS95oHEY9WNOFrZYkynSRQ+48ShSSRe7ZO+mmVM1tEV0B912SQJdGadZv16yyt2cV9li03sp4lEMUsk2h8chX7GLJGQ3zVqMxIO9P5Hoz/6HU0ktL+994RtaP5REpf36SYS4/O3K86TSCSScHwP26O27BhqDrdfsE3phv9luotek0gYhmG6E1kyThbdkEe7DU4kDMMwYUCRJCgGV75GfT0NTiQMwzBhwCpLsBlcdlhj6JKEEwnDMEwYUExubfEcSQ9k3MShGHSqZZK8/DuxGtypk2I1uaYGUZUlk0tQOnkJAO5G8TPumu/Edp1YurWp9oTQ9jaKxYY6MylIJzQVW7zYtoqT3i3LiAorWx+xcqHNkSS0LfGJQttq167TZhdPK1u8xbCfTqZbrOJ+OJPsmm2MSBer/I0YIMaVZBeVYhSNGgpAE5n4psv4SLuxWTxGp5q8hv0A4CbbaCbnEv3MyQZR2FHXoFUMaiaUSZwet7hOeq421mnFI3Um8lVnP/G8CEZpR/F5jUUCikVcZ5xDPNfiE3XOZyK6oMICu008t5Id4nnSL0G7zs6iwOTWFmInk0TVcyQ7d+7E5Zdfjry8PMycORPbtm0DAJSUlGD27NmYMWMGZs+ejdLS0sgGyjAMY4KElh/Y9l6xk0ai6IpEVVXceeedKCoqwplnnon9+/fj6quvxkUXXYR7770Xc+fORV5eHjZt2oSlS5fi+eefj3TIDMMw7dKbJtuj6opElmXU1dUBAOrq6pCWlobq6mrs27cPubm5AFrKze7btw9VVVVGq2IYhokorXMkRq9YIWquSCRJwmOPPYaFCxfC4XCgvr4ef/nLX+ByuZCeng5Fabm3qSgK0tLS4HK5kJKSYrJWhmGYyNCbVFtRc0Xi9Xrxl7/8BWvWrMHOnTuxdu1a3HHHHWhoaIh0aAzDMB1GPn1rq70XFfH0ZKLmiuTLL79ERUUFxo8fDwAYP3484uPjYbfbcfz4cfh8PiiKAp/Ph4qKCmRkZEQ4YoZhmPZh+W8EGDBgAI4dO4ZvvvkGw4YNw6FDh3DixAkMHToU2dnZKC4uRl5eHoqLi5Gdnd3h21qjByWj1t0iz3T1dQh9DiIJTIwTh8VF5MEf7NUasTW7RelnA5G4xvdNE9rU34jKfZvrtHNAVELsaagV2l53o9CmnmC+ZrEfACx+Y18sC5EQN1WTfe87QLNORx9RrutMEddB5apU7msjxyNv3EDNNoYkibLlpDjjU9lKTR51/hu0km82lf8q5FYEvTXhIfsl6+hyqMTYdcottKuJNLfBIx6fJq9WUkwlwzU6EuG2VNaL50VFbZPOUv2FVkYykZKT8TtFzv9mEmflKa3E2EfGy0v3tV7cj/pacaysdq1/14wx4j+YVM5LjxGVX/exd91PYm+abI+aRNK/f38sW7YMt99+e8BQ78EHH0RycjKWLVuGgoICrFmzBk6nE4WFhRGOlmEYxhgZJl5bnVin2+3GAw88gPfffx92ux1jx47FihUrOhtilxE1iQQAZs2ahVmzZmneHz58OF555ZUIRMQwDNM5bDIMJ9ttncgkDz30EOx2O7Zu3QpJknDixAnzD3UDUZVIGIZhYoVg3X9dLhd8pJSF0+mE0+kU3quvr8fGjRvxzjvvBO7a9OvXr0tj7iycSBiGYcJAsHMk8+bNQ1mZWLckPz8fixYtEt47cuQIkpOT8cQTT2D37t1ISEjA7bffjgkTJnR98B2EEwnDMEwYCPaKpKioSPeKhOL1enHkyBGcffbZuOuuu/Dpp5/illtuwfbt25GYmKhZvjvhRMIwDBMGgr0iCfZRhszMTFgsloDLx5gxY9C3b1+UlJRg9OjRoQccAr0ykaQ5RXlqGpGr1hEpI5UHjzlLlEYCwKkmUapI5Y5uss6TxIE4MUmUV7qbtPc+aXnZY6WVYgzHS4S2j8iB9RyFPU1iHGYOwn0HDhLaicmiDBfQuv3GxYsOq3T8qbw6noz3Pw6ITsoAcM5A0X12IBm/DHJM+8aJ64yzaGc665rFY+QkUlC6jr52cR2SV5SnehStk2yNW5SbpidYDPupPJiemwDgJhJWt1Nse0j/QCLHbvaKTsp6NHuNXYpTiRNvf6fWsZlyjBz3L12inL2GyMKTSNw2nWNIzyUaJ90PSlO8sWt0R+jq50hSUlIwceJE/POf/8QFF1yAkpISVFZWYujQoaEF2gX0ykTCMAwTbiySrHl+ifZ3lPvuuw933303CgsLYbFY8Mc//lH3Nlh3w4mEYRgmDEiKBMngssOorz0GDx6MF154IZSwwgInEoZhmDAgKxJkg2Rh1NfT4ETCMAwTBiRZgqS0f/uKVgbtyXAiYRiGCQcmt7ZiybWREwnDMEwYkGWTW1t8RdLzGJgUj76n3UVrm4wdVm3kcnREuviwT5ZPdBAFgH3lonQx3iYOLZUHUymvhUgZ0/uLDsWA1tU1LkGU3jYRqS6V/waDNV7c1zPGDxPa+T89Q2j3c2glrodrRAnmSx8eFuMiDqzHD9cIbSUrWWh/d0Jbk2YQcXA2k3AnESkvdeEFtA7BxPwX9R7xM2Q3YCNy3wSdZwiqmsS43vpalDY740T5KZUpu3Sceul4niLboHJq6m6dqnMMqbOxncRBXXSpa3ENaVPnZL33sjNE9RGVEH9TIUrVJ2b11ayTfjUdVjHuRPK9pDHoycI7i2JVoFi1DsVt+2OFXpNIGIZhupOWORID1RZfkTAMwzBGtMh/DSbbeY6EYRiGMUIykf9yIgkT7RVtKSkpQUFBAU6ePInk5GQUFhYiKysr0uEyDMO0iyRJhrevJK6QGB7aK9py7733Yu7cucjLy8OmTZuwdOlSPP/88xGOlmEYpn0UmwLFZjDZbtDX04iaRNJe0ZbKykrs27cP69evBwDk5uZixYoVqKqq6lDd9tKqetScVrMcOH5K6KP1pZOIiiWeqCv66NQHp6oW4rUHr0l96ioS02GdutseokRy14lKMWqw2CdTVFg5+2lVLldfcqbYJjWvm320xrXYpnXNAWBMvGgCmHjBD4T2v8tElRZVYH1bKapzqOoIAI5Wi0qu6SNFI80UYr5Ha7Tr/TMYT46hT2ff2kLrfVc2im0959cEoiIaT8wnX/l3udA+d5DYr2dUWEUMQmlN9sZqquISz189I0P6nejfR1QIDiIGiuMyxTipau5wjVZBmEnqwCeR75WdzC+kxIvfSz21mcNECUUVmlTxRtuhwJPtEaC9oi1xcXFIT0+HorScIIqiIC0tDS6Xq0OJhGEYpjtpmSPpHZPtXSeaDpG2RVv+9re/YfHixVi0aBEaGrTPEDAMw0Q7raaNRq9YIWquSNor2hIXF4fjx4/D5/NBURT4fD5UVFQEXQyGYRgmEsiSZPj0Or3d2pOJmiuStkVbAASKtmRlZSE7OxvFxcUAgOLiYmRnZ/NtLYZhohpJkU1fsULUXJEA7RdtWbZsGQoKCrBmzRo4nU4UFhZGOlSGYRhDFKsMxdZ+slCsnEjCQntFW4YPH45XXnklAhExDMN0EpMn29n9tweS7LAFjBLPHSxKFakZ3MdfiUZ6l5wn1imnpniAVo56ksh3G8lnmt2iDNFCJK4JOjLPBFJSs1/fgUL7/BFinfcZZ4iS2LQE7eH+l0uUHQ+0iXEfbhJltG4iC9UzP6TfHWqMV0fGgtZwp3x48ITmPSrT/PDISaE9LDVBaFMzPr37031MdP1nEyNNKg+mhoFUKg0Ap5rF8RqTLsbpGS3O/W3dXyG0HTox7impEtonickllZrTcy99iPh9AICh/cS46DFKTRClt1TuayUnwZAkUeoLAOmkzjs1zbSSH1rar/cbTRTZqHGL5zOtb0/NJhO68NkOdv9lGIZhQkKSTZ5s50TCMAzDGCErsuFzJEZ9PQ1OJAzDMGGgNz3ZHjspkWEYJoqQrApkq6Xdl9TJwlZPPPEERo4cia+//rqLI+48fEXCMAwTBmTZ5NaW3PH/47/44gv8+9//RmZmZiihdTmcSBiGYcKA2UOHrX0ulws+n6ikczqdcBKVZnNzM5YvX46HH34Y1113XdcHHAK9JpEMctrR6GmRslIp4ohzxLY0WawZ7o8n8kgdV1jVIsojm4j0c3+lWGvbQ/p/kCx+vl+89rJXbqwWt2kV5ag0BsknusCWNYknKwAMTRZdXU/4xO2+/Y0ohR4zQHT21YNKgj1EcjmGONpSKe4334ly7P8YlqrZRiqRjvYlTrBuHVlyW2gNd0ArN7Ur4ljQeusDE0l9dfL5ap3x7kscbulnzkwVz0U/cTV+8aMjmnVmEBfdKjJ+FnIL5YzhoivEsP6JmnVOGCJ+B84kcuAku7hOKnV21Ynnnh7JpD69kzy8F08e2LOQOQW96Qf6zax1i3GmOsRtnnKL54mtCyfAg00k8+bNQ1lZmdCXn5+PRYsWCe/9+c9/xqxZszB48OAui7Gr6DWJhGEYpjtpKWxlkEhO/wNVVFSke0XSlk8++QR79+7F4sWLuz7QLoATCcMwTBiQbVbINqthP4CgDGj37NmDb775Bj/96U8BAMeOHcOCBQvw4IMP4oILLuiagEOAEwnDMEwYkCTJcEK9I6V2b7rpJtx0002B9rRp0/Dkk0/izDPPNPhU98GJhGEYJgwEO0cSC3AiYRiGCQPhTCQ7duzo9GfDAScShmGYMCBJsslkO1+R9Dh+nJkA9bRsV6mvFDuJVFTyiY6gUlOd2PZr3X8lT6PQdiSIktWxaaJsVm6qIZ8n26gUJZwAoNpECSYUUfIqn6oV2r5EUTpqkbVyVA9x0W3wiGMxLkNUj1DZst4/VTVEUkm3MdApSo5ddW6hff5wcewSdRxZE23iqesgElcq5VXI/Wi7jrtyHyI/pRLXRLfosuuVRLflhIbjYlvnh0L1ivvuVcTxHWgVx0LuLx7zS84ZoFlnSZXo9jsiXZTzfkackc8ix3TW2dp1jk4Q1yk3HBXa3vhhQrveKx5jVRXPTXoOANrxTbCSY0bkvnRGwQKtxFvyiuPXXxaX6Ufkv74+4vHoyLyFGZKJjTyX2mUYhmEMka0WwNb+T6xsjZ2f36i8tqJeMiUlJZg9ezZmzJiB2bNno7S0NLIBMgzDmCDJsukrVoi6PdHzkrn33nsxd+5cbN26FXPnzsXSpUsjGCHDMIw5rbe22n/Fzq2tqEokrV4y9957b+BeZWVlJfbt24fc3FwAQG5uLvbt24eqqiqjVTEMw0QUSTZKInxFEjb0vGRcLhfS09OhnPY9UhQFaWlpcLlckQqTYRjGFL61FQFavWTmzp0b6VAYhmFCRlYU01esEDWygfa8ZJYsWYLjx4/D5/NBURT4fD5UVFQE5U/TFnn/P4DmFkmjZBMlf7JTdEP1nxQdb+UBotRRI8Nt5z2hn8gK5UZR/qtaRQdXmUiOAUD1iNLGmvg0od1HPSW0LSdFyWYGcQcGAL9DdOKVmkQJ8YD+ohT3RIMofa5q1EqKAXGZjERxu33jRQlmZh+xP84ijhV1fQUAm4m816dSOar4eb3b09R91la+V1yHTL4uVvE8Um2iG7PcILo1A4BcT27JJoouuVR6nimLMtyfDRO3AQBVw/pq3mvLU0R6O3esWMtiWO2Xms/4q0TJsJQift/kZlGenkikzj67eEyDOYaa/mYiiSffCb3vCBpI3BbxXKPfM4tMfswVKzDobMO4gkWyKpANVFudLWwVjUTNFclNN92Ef/zjH9ixYwd27NiBAQMGYN26dbj00kuRnZ2N4uJiAEBxcTGys7ORkpJiskaGYZjI0ZtubUXNFYkRy5YtQ0FBAdasWQOn04nCwsJIh8QwDGOIJMv6T+y27Y8RojaRtPWSGT58OF555ZUIRsMwDNMxJEUyTiQxJP+N2kTCMAzTk5FkGTDy2uIrEoZhGMYISbECFptxf4zQaxKJv/YE1MYWVZPaLNZPpyouJVU0sVOowqpZNGgEANVG6r4TE0eZGkUS5FMnxPVZ4jTLqFZRCZN8/DOxn5g40piUWq2KyB8nGvjRE5/46KFvHDFH1FHj9HOIp1WSXfzPi9Ypp7XlVaKkkVSd+usqUYtRg0S/2C95RPWT5WS5ZpX+GvEYIEk0ZQRVAJEfAj9RxUlurfEmLKKiTWMgStVmdlENSE0JASDFKp7P1BC04DzRvFM+9K7Q9tSJSicAsA4eIcZB+qVmcTxVu2gUmSTROHX++yZmqVTJKDeLKkT/sRJxm1RxBcDfJI65RM5nf534HfA3ifshxfeB0kWqLZhckRj29TB6TSJhGIbpTiRFAQyeFZH4ORKGYRjGEEkGdK6ahP4YgRMJwzBMOJAV40Ri1NfD4ETCMAwTBli1xTAMw4SGiWoLrNpiGIZhDDGzQeErkp6H6vVC9bYY4lEjN0v/gULbT+TBvspjQluK0xrn+U+JEkp/jSjrVMg2qORYSiAyXB0kr3FcmhPTL8or1T7JmnVavv1Y/Ei9aIRnHfADcXliTJig81+V7CZmeqfEOKgxoZ9IXKmMWW9Skspmm/fvEWNIEs0mKW46dgBkRx+h7auuELd55ngxLGIaaFFF+arkF6W+ACDVioag9FykqDViDJKzn2YZhZwX/vJDQrvijceEdup5Y8UV+LXGm55vxfcsA4aI2yQSeFnHEFRcQGJchg4AACAASURBVE/+K54Haq0ov6bfQ4nMKdDjA2jlvPCK0nJ6hUDHX7J04U+ibDLZ3sFEUl1djTvvvBOHDx+GzWbD0KFDsXz58qjwHYydlMgwDBNNtE62G706gCRJuOGGG7B161a8/vrrGDx4MB5++OEwBd8xOJEwDMOEAUlWICkGrw4mkuTkZEycODHQHjt2LMrLtQ/WRoJec2uLYRimW5ElkyfbWxweXC4XfD7xVqLT6YTT2f7tbr/fj5deegnTpk3rklBDhRMJwzBMGJAUm8YSR+D0XOC8efNQVlYmdOXn52PRokXtfnTFihVwOBy45ppruiTWUOFEwjAMEw6CvCIpKirSvSJpj8LCQnz77bd48sknIUeJ8itqEomRIqGkpAQFBQU4efIkkpOTUVhYiKysrEiHzDAM0y5SkE+2d6Rs+KOPPorPP/8cTz31FGw2g2dUuhlJVWk168hw8uRJfPXVV4HJpMLCQtTU1OCBBx7Atddei6uuugp5eXnYtGkTXn31VTz//PMdWn/T9vVQG0/LNYnc0ZIhSlw1skJyMsg6Mlo/dVAlskMq/1WJ9JF+Xk9/LsWLMlnfd+LlsEbaSCTGfh25JJXJaurZm8iSVR3pqOpu0lmyTX+TjituRyHbpeNJpby0n+4nAHhcpUI7fsz54gIJpDY6cfeViDuzv5bUZ4fWndaSNlhcgEidmw9+KrRlcg60vCl+xnP4a6Fdc0g8Txxp4n7YM8VzEwBUtyjvlazGP1pyovidoOcqlfrqoZL/yr2uknaWPL0NvR9pkwlsyU5k91T+G98H9qldc7vI+8U7gI5TeABbPCznTAl6fQcOHEBubi6ysrIQF9eyH4MGDcLq1atDDTVkouaKRE+R8NJLL6GyshL79u3D+vXrAQC5ublYsWIFqqqqokI/zTAMo4tkcmtL6liFxBEjRuCrr74KMajwEDWJpC1tFQkulwvp6elQTlsuK4qCtLQ0uFwuTiQMw0QtLVdxOrV0WjG5yutJRMdMDSHaFAkMwzAdpvXJ9nZfUfnz2ymi7oqEKhIyMjJw/Phx+Hw+KIoCn8+HioqKDk1QMQzDdDeSZFIhMYbqkUTVnrQqElavXh1QJKSmpiI7OxvFxcUAgOLiYmRnZ/NtLYZhopsutkiJZqLmiuTAgQN48sknkZWVhTlz5gD4XpGwbNkyFBQUYM2aNXA6nSgsLIxwtAzDMCZIkvFVRwcn26OZqJH/hpvalx+GetqhV/WJE2CyVcynzXWiRNORrnVcpfjcbqGt2EUpKHX3lYmDsPuIKHX01Gtlg44holRUI7Mll9HUCVXWcZqlMk0qhzSTzdJ+PaiUVCMFNfvPLIh7yab1r4k0mo4/AMh9RFmsSiTcVOLqqxadfKkk2d9AXJB1ljE7hhrZrY7cmo4fXaevUTwPJEXchq+JOOQCsGcQuTp14nUQWTgZK5U6T9NzANC4U0tW8dxzVxCnZBI3/R7rLdNcK36X6WdsTvH8lxOT4Zx7pzbWTuA//JnWfbgtFhvkIed2ybYiTdRckTAMw8QSqmwBZAPVlhw7P7+xsycMwzDRRgzdvjKCEwnDMEw4kGVA5QqJDMMwTCdRJdlksp0TCcMwDGMEJ5LYo7mmHv7aFhUNVXbQtqdeVKhU7y8V2tYEreEfVakkDk4XF6ivFZcn6ihrqqgMszi1ag8zM0Sq8FGIkoma4umukyqAPESNYxIDAK0ySdOmE5DGyi9vvXab1j6i6spLlEdKHDGwlImSiRwPQEeFRW89+EtJW9wvqoDTGyuqTKJtzTGiSrBGreGlxsiRxG1xJonb8Ijjrafa8pw4Lm7XI9bVsJI4PSdF01FbuvjAsJ7pIx0f2rYlicab7qoaoS3btD9fFgcxKm02qAcCoKGiWmgrbgnGNqUdQFYAI1EsP0fCMAzDGMG3thiGYZgQMUkk0WUsEhKcSBiGYcKBJBnLf2NIGsyJhGEYJhz0IosUTiQMwzBhoGWOxGCynedIGIZhGEMUi/FVB6u2eh4JGf2hJrVIEKnBIjVyU4gks7FSlB02VWqlo44M0daeyiUVIhn268hPBfTMEMmJp5GOEqkuNcXTRWMiSGp1Ewmxn5r3WXSqvNEviJ7RoAFmppp6aCTcdaJRIZUD60FlxmYycYqFyoF1xsZTLcpNLYmJ4mesWmPNtih907Rv0vE1Oe40LptdR85OpLiK3VjuKxHJsWp2fgPa84TWTyfnnj2NjE0Q57e9Hyk3QeK0p4jSaMkhtkOCb20xDMMwoSEDhrkidm5tBb0nb775pu77b731VpcFwzAMEyuokmz6ihWC3pN77rlH9/2lS5d2WTBGlJSUYPbs2ZgxYwZmz56N0tLSbtkuwzBMp2h9INHo1UGi9XfQdE+OHDmCI0eOQFXVwN+tr/feey9QEjfc3HvvvZg7dy62bt2KuXPndlsCYxiG6RQSvn+WRPfV8VVG6++gaSKZPn06Lr74YjQ2NmL69OnC684778SiRYvCHmRlZSX27duH3NxcAEBubi727duHqqqqsG+bYRimU8gWqAavjha2Csfv4Lp163TfX79+fYfWY7on+/fvBwBcc8012LBhQ4dW3lW4XC6kp6cHTAgVRUFaWhpcLhdSUlJMPs0wDBMBJBmA0XMkLZckLpcLPmKC6XQ64XSK9pHh+B1cvXo1FixYoHl/7dq1+NWvfhX0eoJOiZFKIl2Ft6kJakOLtJW6hirxopOspIjSUUcQtaJpnXclQZR1ampeU2kolXAGoTHXuMvSmuydWCeVXFKnWCmIdZjVCKdQh1atO7AWbwORp1K3X5O65BYdB+emk2KNdcVEdpwwsL/QlhPEL75EjweAuMwsoU3HV20iNca9xv2AnmOwcS10urxpvXsAfi89D8Txpd8pKrHXPweMXZ/p90zjYqxznsgJomOwlzg6K6Sfjm9XokqSURoJJJJ58+ahrKxM6MrPzw/r3Z73338fAOD3+/HBBx9AbeNSfPToUSQkJLT3UV2CTiRerxcvvvgi9uzZg+rqamHDRUVFHdpoR8nIyMDx48fh8/mgKAp8Ph8qKiqQkZFh/mGGYZgIoKrGLvKtFBUV6V6RULryd7BVPOV2u3H33XcH3pckCf3798fvfve7Dq0v6ETy4IMP4oMPPsAvfvELPPbYY/jNb36Dl156CZdddlmHNtgZUlNTkZ2djeLiYuTl5aG4uBjZ2dl8W4thmKjFr6rCP9yU1rn2YBNBV/4O7tixAwBw55134o9//GOHP08JWn+2bds2PP3007juuuugKAquu+46rF69Grt37w45iGBYtmwZNmzYgBkzZmDDhg247777umW7DMMwncGnmr86Slf/DnZFEgE6cEXS1NQUyJxxcXFobGzE8OHDsW/fvi4JxIzhw4fjlVde6ZZtMQzDhEyQVyQdIVp/B4NOJMOHD8fevXtx7rnnYtSoUXj88ceRmJiI9PR08w8zDMP0Mvxqy6s3EPStrbvvvhsWS0veKSgowL59+/D3v/8dK1asCFtwDMMwPRU1iFesYHhF0ioR03vv5ptvBgB4POGTz3UlDccq4atpeXDHQqSiVDpqTYgX2tR91ucjLrsA+gzNFNoaua+pc297kRtApbhUykjkp/4mUZYL6MghqTTUJE7qBgxo5dSqVxyv5pOnxG0QaWgwTr0UKu81g0qOAaDPEPHqmspPqaMzdRiWiCTZPnCINk4iRwUZG43cmoy/ZNNx6q0TnXhV6uhMttFE3Kwt8XbNOqm81+IQvxN0vKlzcjDHkB4DOr7eJmMJse4xrxVdh+l+NFUeNYxBdjZD3NPO05uuSAwTSXv+Wm2RJAn/93//12UBMQzDxAKqyRxJLOUYw0TSKhFjGIZhOoYfJsqs2ClHwvVIGIZhwoHZrS2jKrw9DU4kDMMwYcD01hYnEoZhGMYI/+lXe8TQnS1OJAzDMOHAzGuLr0h6ILLVApx2KPUQqSJ1LqXySIrNqXXG1EguiTMpbdMHeOjng3HZpfJeyU6koWQdit1c2CjHEclwTaXQbvxOrH1ApdIAUPX5l0I7LlU0oPMRmWfSGUOFtveUKA+mx0cPi1XcdyolleLMHZw1DsEOcZ1UBk6luNRF118vSlFb3iTbJfJUegz15L4UibgOK4nJQlttEONwkOU1LtIApHjxHFcbibs1kffSsaNtPfdfKs2lxyRhyCDxA+R8tqRqx8ZddlhcJxlvxap1ZBY2YQ3iexckqqrCb5AtZE4kDMMwjBFmflo82c4wDMMYwre2GIZhmJDwQ4Xf4LFD8/JtPQdOJAzDMOHArLAVX5EwDMMwRvhh/EAiX5F0Mffddx/ef/992Gw2OBwO3HPPPRg9ejQAoLGxEUuWLMEXX3wBRVFw1113YerUqR3eRkJmP6hJLYoNquRwV4u1uqv2fyu0nVliBbOaQ2J9ZQAYNON8oU3N96gSiZolgiyvZ4Yom6iEqLKGKoCUVPNKbO6S/UK77vBxsd/EcBEAvPWNQruRqOBk8pmEAalC20LLjAZRw13uIyqVLGSsNLXQdcbXVy+eB7QuvFWn/GlbNCoui9a4UFNznaj1NEacFJ1+JYnERdR8tI65TBRZUpxWheje/7HQ9hAzRKraqneJ6j5qgKmnhKQqONmkdrxM1Gh6YxF/1rlC208MLS39BxpuA7ausmwEfKoKH6u2uo/Jkyfj7rvvhtVqxc6dO3HHHXfg7bffBgCsW7cOCQkJ2L59O0pLSzFv3jxs27atw8XpGYZhupPeNNkedD2ScDJ16lRYT+u7x44di2PHjsF/+r/QN998E3PmzAEAZGVlYdSoUXj33XcjFivDMEwwqGpL3fb2XrGUSKLiiqQtRUVFuPDCCyGfflipvLwcAwd+fzmakZGBY8eORSo8hmGYoPD5W17t9seQR0q3JJIrrrgC5eXlun3vvfcelNP3Rrds2YLXX38dRUVF3REWwzBM2PCbPNlu1NfT6JZE8tprr5kus337djz66KN47rnn0K9fv8D7mZmZKCsrQ0pKCgDA5XJh4sSJYYuVYRimK/CbTLbHUiKJijmSnTt34sEHH8S6deswaJDor5OTk4OXX34ZAFBaWoq9e/di0qRJkQiTYRgmaDx+FR6fwauL6/Ded999yMnJwaxZszBnzhzs3bu3S9dvRFTMkSxZsgRWqxW33XZb4L3nnnsOffv2xYIFC1BQUIDp06dDlmUsX74ciYmJHd6G0qcvVEvL7npJ3ey41CShTeWotaUuod1/7AjtBogU0VMryiGtiaLKTE5OE9q+SnEbVOoLaOt3a0LoK66TmuLpfsaZIrQVu1i/O56MTX2ZOHaeerGuNgA014uGfQlpovzUQowem+tESSytD27J/IFmG2ZQw0SV1LOn8mxAf8zbQsefylE1Roc6ZohUmqv6iYSVmnfaiExcTx5MjR+pmSeRIVOZsmrVSl6t4y8W39hdLDQ91dXiOomkm7bj0/pqtkHrpXuIbFwzfnRsrDp14alRaeoA7TJtt0Hl17QdAt19a8tI/RpuoiKRfPDBB+32ORwOrFq1qhujYRiGCZ3uvrXV9vm6tupXOYh/KEMlKhIJwzBMrGFWare1z+VywecTr4ScTiecJg/AGkHVr+GGEwnDMEwY8PlV+AwySWvfvHnzUFYmumXk5+dj0aJFwnvRrH7lRMIwDBMGVBjPkainXRuLiop0r0gooahfww0nEoZhmDDg8cNQmWU5/bBiRoa5B14wtKpf169fr1G/hhtOJAzDMGHAb3Jry9/F8l8j9Wu46TWJROmXCTS3yAsV4gDqI3LgfhNEGafy+ddCO/GsszXr9zc1CG37YFGyWv+16KrroLW5rcYSTQAaaSOVeVoGDBH77Q4xRofOCeUXJZj2MaIUVDl6UGgPGjxYaFf/+3PtKkm9dEpzrSj31dT7douSYolIowHAPurHQttz9JDQtmRkiTGdEl1gqStsy5tUXkpktNTdl9Q+l4mLrqrj/uunx9BLapvT+utE0q06xTYASFVHxWXIbRLqEg2FSJAV858B69Bssg3xuNuIdFeOE889OUF7q4Y6MFuctM57x+unU8mwSmTgIMeEOiN3Jd0t/zVSv4abXpNIGIZhuhOzmu1GfT0NTiQMwzBhQDW5IlFjyCKFEwnDMEwYaLVIabdf4UTCMAzDGMDuvwzDMExI9Cb3X04kDMMwYcDvVw0lvl0t/40kvSaRSKmZgO+046sk+s8og88RF/5WtF/ulyo+MCQ7dCSDRE4qEXmv87zzhbbnsCgplhPE5fVcS2XiLosEUc7rS84UY3CLMlu63wCgUklwguh8jEHnCk1Lg+j62n/iTM06pcPi+PnI2Pi+E+0gGo+L8muVlJXzN5L9AOCrFKtkWoeMFNdB9tVC5Kd+4twLAJ5vRYk2lVtr3ZVJP3FSVj1aZ2S5D3EMthg7OsMiujGrqrbknpoqyr7pvmt+rqj816qVmqsKkckOEuOwNInHRKmsEtrWQWeIn08i5xUAxInfI7WaVD6lPlHUwdmrdXD21xDX7RHjtNttuw6ruF90bEKBVVsMwzBMSKgwvn0VQ3mEEwnDMEw4aPb50WxQtN2or6cRFRUSW9m9ezeys7OxYcOGwHuNjY34zW9+g+nTpyMnJwc7d+6MYIQMwzDB0WqR0t6L50jCwKlTp/Dwww9j8uTJwvvr1q1DQkICtm/fjtLSUsybNw/btm1DQkJCO2tiGIaJPMHayMcCUXNFsnLlSixYsEBjMPbmm29izpw5AICsrCyMGjUK7777biRCZBiGCRq/anJFEkPy36hIJO+88w5qa2uRk5Oj6SsvL8fAgd+bLGZkZODYsWOa5RiGYaIJoyRidrXS0+iWW1tGlb3eeustPPLII1i/fn1YY/Ck/iDgbVNRL7rTWmRJaCeNFmW09loxdt3DT2Szco0ocZWbRXdUjbyXOM+q8TpuqVS+K4uHr8EqfqZZFuWVdov2/wY7XWWzKOuUSFu1JxrHBABnilJnjYyWfCa54oB2HW3wHdmveU8eILore1KzhLbkFaW3SrXokIv+QzXrtBAHZz9xjqVyX41DM5HmSlRaCkDVyHnJ2SSJ56JqIw7OVtGdGQDkphrxDZt429fTRzzXFJ84NvWqVvIabxWPkc8i7qstS5SF9yHSaH9fsR6G26mtuUH/I1cGnCW0rdVHhLbUVCe2Za07sNKf1OEgx8TXh7gpU+lzF5am9XhVNHvbn1D3eDmRdAijyl4fffQRvvvuO/z85z8HAFRXV2Pnzp04efIk8vPzkZmZibKyMqSktGj0XS4XJk6c2B1hMwzDdBqfajJHEkO3tiI+2T5hwgS8//77gXZBQQFGjRqFa665BgCQk5ODl19+GaNHj0ZpaSn27t2LRx55JFLhMgzDBEVvmmyPeCIxY8GCBSgoKMD06dMhyzKWL1+OxMRE8w8yDMNEkO6ukBhJoi6RrFy5Umg7HA6sWrUqQtEwDMN0Dq9fhdcgWRj19TSiLpEwDMPEAn6TOZJYkv/2mkTyTbU7YEnQRJQUVqLUqGsW2/HWdKHdx6ZVdlADthq7aKQ3IEUc6vhKsRa6SlQxqqxzaIiip94hKlBq3OJ+nWgQ1WlJcVqVi5Uo1hxWUSUUR+qQW8XFNeooAGhWRGUS/TJJEFfiTxfVOopXrLstJWkNLD12UZHmIUp2xU7UUWkjxBX4tIZ/ykCi5qv9TrNMyJBzzW8jKjg7qftOzgOJjE3Lm+I6qdLO0igaKqpk7BwKOajQKhPpd6YxcaDQ7kPMP+tk8TyqqdcxWDT5HU1KELcRnyTGafWKSkgA8FhEVdspjxg3/ebWNov9FllClnFYQdObLFJ6TSJhGIbpTvx+4wl1f+zkkeh4IJFhGCbW8Jk82R4u+a+eZ2G44SsShmGYMBAJ1VZ7noXhhq9IGIZhwoDP7zd9dTXteRaGG74iYRiGCQPBPpDocrng84kWSU6nE06n1ibJiLaehX//+987HG8ocCJhGIYJA80+P9wGXlutqq158+ahrEz05svPz8eiRYuE96LBs7A9ek0iOdnkCUgYtf8kiP8NnGoWZYYpDtHUrlmn2HK9R1xHTZMoJXUTg7ZBKWJNa0f9caEt6chTfcQYsonEUdssxhBHTBr1JvdU8SPwq3QsxH66To9fewop5MtjI/JSWnfcSyTIFlmU7p6C1lTQSobHIpN9o0pn2WrcBnAqSTSCVJLFtkTk1/RHgp4Xej8hZrYYso+MhSq243Tqq9vjxNsYNg8x3qRmh17xoPr6iPJ2QLsvbrpvZDcsxCiy3k3PI+1oKGQ86flJY4i3iMt7/dpjaFPE7drJZ6qaxH4qa7YpXXe3P1jVVlFRke4VCSUUz8Jw02sSCcMwTHcS7K2tjAytM3JHMfMsDDecSBiGYcIAu/8yDMMwIRFJ00bqWRhuOJEwDMOEAa/XD6/BZLtRX0+DEwnDMEwY8Kuq4VUHmzYyDMMwxqiqtpQy6Y8VoiaRvPDCCygqKoLVaoWiKNi4cSMAoLGxEUuWLMEXX3wBRVFw1113YerUqR1ef1tNN1XvUvNTuyIOC73PWeslmlkADUT+ayUr9Zg8xeqPE+V+kqpdvkYlslgdSWVb6GReDXED1iPVoZVUtkUmg0fllYD2+2GkpQcAhSzfSBxbFVm7DSoN1VFkk42ITb2x85Dj3EjiPkXclc2OqR4eEiiVU8eRWunUnbnBox0Lh+YzogOu4hCdeOkaHNBSR+S7dHwtmrjEsTh+StRnVzXquC3rHNe2aJ2pxYNIxwoAPH7xM1S+XtcsfgeqSVzxVhlnG0YVPKofUA2uSHS+4j2WqEgk27Ztw1tvvYW//vWvSExMxHfffW/fvW7dOiQkJGD79u0oLS3FvHnzsG3bNiQkJBiskWEYJrL4/Sa3tmKosFVUeG09++yzyM/PD5TQ7d+/f6DvzTffxJw5cwAAWVlZGDVqFN59992IxMkwDBMsLVckxq9YISoSyaFDh/Dpp59izpw5uPLKK/G///u/gb7y8nIMHPh9gZuMjAwcO3YsEmEyDMMEjc/vh89n8IqhgiTdcmvLyCPmvffeg8/ng8vlwosvvojq6mpcffXV+MEPfoDzzjuvO8JjGIbpevyq4RyJaYnIHkS3JBIjjxgAyMzMRG5uLmRZRmpqKn7yk5/gs88+w3nnnYfMzEyUlZUhJSUFQItT5sSJE7sjbIZhmE6jmiQSwyTTw4iKW1u5ubnYtWsXAKChoQEff/wxzjqrpY53Tk4OXn75ZQBAaWkp9u7di0mTJkUsVoZhmGDwq6rpK1aICtXWL3/5S/z+97/HZZddBgDIy8vD+eefDwBYsGABCgoKMH36dMiyjOXLlwcm5TuCx/+9vDPRJsoIqQup2yfeu4zzUcdbHSdTk5RMVgkXkUc2+8SYbIreoRElmTRuCpWaUvkkACSQ9/rYxR2xEQlmM/kviroaA1oX4nrSplJpKh2l/CA5XvOehexKX7v4ximyzoYGsf1dvVaOeopIQ5PixGNgV8Rt0PGlkm8q7QUARSLyXyJhlYk4lx5jKlEGtDJlYl4NmayD/oC5fXoyWmP3X+roTOXUVO6rJ/VNT7Rp3mtLea1baNPzRk/yTfeNSvfpd5uOjamMvAOoqskVCSeSriUuLg4PPfSQbp/D4cCqVau6OSKGYZjQ4OdIGIZhmJDwn1ZnGfXHCpxIGIZhwoCqGl91xNCdLU4kDMMw4aA3PdnOiYRhGCYM9Cb5LycShmGYMMCJJAYZkGBH8+nJLSrVpTJOCpUIWmXzx2+orPM4sSE90SBuk8ojqaMroJUQJ8eLh6+ZSHGpTLlCR/JKoZJXKhmmMVDJLKB1Hf6mqkFoU8dVuu/J8aIDMXWB1aOCOsMS6S2Vjta4deImX2xPgxinwyrufN84Y6dkPewWY7kvHbsmcm7S5QHAS12LSdtMJt7s17pZU4NrOn703KLnRR8ix27ScYCm8l4KPS8GJcaJ29SZZKhpEo8rPc5W8uXvR9yu7WY6/g7g95tMtrNFCsMwDGMEy38ZhmGYkODJdoZhGCYkVJMKifxkO8MwDGNIJCxS2qs0G244kTAMw4SB7r61ZVRpNtz0mkRS2egO1AKnSg6q2nLr1GRvi6yjIqLKGKrOSbKLQ60xiyMnFVWf6H2GKmfM1GdUsQIAaQmicR41vaPKmX0VdSQG7TY9RKlCt5tKtklrutMv2JGaRs02+hJlF1WbaVVGxGBR5xhSZVGt5jwR+6npID0HqLGhHlR5JJGa49r69drx1jlVBNw+4/NdT/1Ej4neudMWaoRKVXPBKB2p6SWNkx4ft0/7PaXHnR5nuo1TxFDUa+m6H3fV64Hf22zQ33HVnxHPPvssbr/9dt1Ks+Gm1yQShmGY7kRVfVB1pNVt+4GWGks+khSdTiecTmeHttdaafbPf/4zmpubMWfOHPziF7/oeOCdgBMJwzBMGFB9fqg6V01t+wFg3rx5KCsrE/ry8/OxaNEi4b1orjTLiYRhGCYMqH6TK5LTfUVFRbpXJJRQKs2GG04kDMMwYSDYRJKRkdEl22utNHveeecFKs1Onz69S9ZtRlSU2i0pKcH8+fORl5eHSy65BI8//nigr7GxEb/5zW8wffp05OTkYOfOnRGMlGEYJjhU1R9IJrqvLn60/Ze//CVcLhcuu+wy/PznP8fMmTMDlWbDTVRckTz00EOYMWMGrrnmGtTX1yM3NxdTpkzBueeei3Xr1iEhIQHbt29HaWkp5s2bh23btiEhISHSYTMMw7SL6m02UW2139cZjCrNhpuoSCSSJKGurkVW2tTUBEmSkJKSAgB48803sXLlSgBAVlYWRo0ahXfffReXXHJJh7ZRVusOSP2oVJRKXqm8l8o6qQwXMJf7amWyZJtknTRGAEi0ie/Fk21SGSf1i9MzzqMyY2rCSD8TT8wR7TrSUZ8qLkPHl0KNIalkM6OPXfMZOhb0kGhrnZNa3TZt/fp+DlHOqzHrREPvygAACItJREFUNJHz0vFv1KlnT+Oi8l4aJ5XhBiO3psvQ/eijs++adZJ19HOI450QxDraoidbpiajGnkvaQdjvEklw/R7Sc/FU2QdVMYcCn6/H36DW1uxZNoYFbe27r77brzxxhuYNGkSpk2bhgULFmDQoEEAgPLycgwcODCwbEZGBo4dOxapUBmGYYKiVf7b/q0t4+fVehLdckViJlt7+eWXkZeXhxtuuAEVFRWYP38+Ro0ahTFjxnRHeAzDMF2O6vebTLbHzhVJtyQSM9naCy+8gLfffhsAkJaWhh/96EfYs2cPxowZg8zMTJSVlQVudblcLkycODHsMTMMw4SC6vOZPEcSO1ckUXFra9CgQdi1axcA4NSpU/j4448xYsQIAEBOTg5efvllAEBpaSn27t2LSZMmRSxWhmGYYPD7WixS2n35zAvN9RSiYrL9wQcfxB/+8Ac8++yz8Hq9uPTSSzFlyhQAwIIFC1BQUIDp06dDlmUsX7484CXDMAwTtZg8RwKjvh5GVCSSUaNG4X/+5390+xwOB1atWtXNETEMw4RGyxxJ+/MgPEfSA4m3KmhVNFKJIJVH0jZVfeq5/1KoVJG66MKkVveJBu1lL32POgZTqNTRxBxY9zOdgbr30u2abcNDxv9Qldb9lx4Ti0LdZql0V+wPwphXV+bdFrP9MPs8oHdeiNBjHMzxoduln6kjjrd6Y0Frl1OpLW1r6t3Ttk7t8o6en3R5PRt2ui9eg5rpAGC3iHJfm6Xr5L+tDyQa9ccKvSaRMAzDdCfBWqTEApxIGIZhwoDf7zN5IJETCcMwDGOA3+uB39O+DYrfy6othmEYxgiTBxLBk+09D+rn1BE0k+06k6gS9WUik6h0gt5svj6Ycs5mNZ87M9mudsVku8l2u2SSmnprycaT7RbFWEChhwTjhVQY74fZ54GOT7abbRPQG2+xTTepNxY09o7GaSHHWK+0cUfPT3pudoUwhG4zlN8JSmZ/p2EiyezfsQqI0YykdsUvB8MwDNNriYon2xmGYZieCycShmEYJiQ4kTAMwzAhwYmEYRiGCQlOJAzDMExIcCJhGIZhQoITCcMwDBMSnEgYhmGYkOBEwjAMw4REr7BIKSkpQUFBAU6ePInk5GQUFhYiKyur2+OYNm0abDYb7HY7AGDx4sWYNGlSt8ZXWFiIrVu3oqysDK+//jrOPPNMAMZjFO742oupvfHqjpiqq6tx55134vDhw7DZbBg6dCiWL1+OlJSUiI2VUUyRHCsAWLhwIY4ePQpZluFwOPD73/8e2dnZET2v2osp0mMVk6i9gPnz56sbN25UVVVVN27cqM6fPz8icUydOlX96quvNO93Z3x79uxRy8vLNbEYxRDu+NqLqb3x6o6Yqqur1Q8++CDQXrlypbpkyRLTbYczLqOYIjlWqqqqtbW1gb+3b9+uXn755abbDndc7cUU6bGKRWI+kZw4cUIdP3686vV6VVVVVa/Xq44fP16trKzs9lj0TuBIxdc2FqMYujO+YBNJJMbsrbfeUq+77rqoGau2MalqdI3Va6+9pl5xxRVRNVatMalqdI1VrBDzt7ZcLhfS09OhKC2unoqiIC0tDS6XCykpKd0ez+LFi6GqKsaPH4/f/va3URGfUQyqqkY0PjpeTqez28fM7/fjpZdewrRp06JmrNrG1Eqkx+qee+7BP//5T6iqimeeeSYqxorG1EqkxyrW4Mn2bqSoqAibN2/Gq6++ClVVsXz58kiHFNVEy3itWLECDocD11xzTUS2rweNKRrG6v7778ff//533HHHHfjjH//Y7dvXQy+maBirWCPmE0lGRgaOHz8On6+lLoDP50NFRQUyMjIiEgsA2Gw2zJ07F//617+iIj6jGCIZn954mcXb1RQWFuLbb7/FY489BlmWo2KsaExAdIxVK5dffjl2796NAQMGRHysaEzV1dVRNVaxQswnktTUVGRnZ6O4uBgAUFxcjOzs7G6/VG1oaEBdXR2AlgI9b7zxBrKzs6MiPqMYIhVfe+NlFm9X8uijj+Lzzz/H6tWrYbPZTLfdHXHpxRTpsaqvr4fL5Qq0d+zYgaSkpIiOVXsx2e32iJ9XsUivKGx16NAhFBQUoLa2Fk6nE4WFhRg2bFi3xnDkyBEsWrQIPp8Pfr8fw4cPx+9+9zukpaV1a3x/+MMfsG3bNpw4cQJ9+/ZFcnIytmzZYhhDuOPTi+nJJ59sd7y6I6YDBw4gNzcXWVlZiIuLAwAMGjQIq1evjthYtRdTQUFBRMfqxIkTWLhwIRobGyHLMpKSknDXXXfhnHPOidhYtReT0+mM6FjFKr0ikTAMwzDhI+ZvbTEMwzDhhRMJwzAMExKcSBiGYZiQ4ETCMAzDhAQnEoZhGCYkOJEwMcu0adPw3nvvRToMhol5OJEwDMMwIcGJhGEYhgkJTiRMzNPc3Iz7778fF1xwAS644ALcf//9aG5uBgDs3r0bkydPxrPPPosf//jHuOCCC/Dqq69GOGKG6VlwImFinrVr1+LTTz/Fpk2bsHnzZuzduxdr1qwJ9J84cQJ1dXV49913cf/992P58uWoqamJYMQM07PgRMLEPK+//jpuvfVWpKamIiUlBbfeeis2b94c6LdYLLj11lthtVoxZcoUOBwOlJSURDBihulZcCJhYp6KigpkZmYG2pmZmaioqAi0k5OTYbF8X+MtPj4eDQ0N3Rojw/RkOJEwMU9aWhrKy8sDbZfLFXB7ZRgmdDiRMDHPZZddhrVr16KqqgpVVVVYvXo1Zs6cGemwGCZmiPma7QyzcOFC1NfXY9asWQCAnJwcLFy4MMJRMUzswPVIGIZhmJDgW1sMwzBMSHAiYRiGYUKCEwnDMAwTEpxIGIZhmJDgRMIwDMOEBCcShmEYJiQ4kTAMwzAhwYmEYRiGCQlOJAzDMExI/D/fjQwDdvvblgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_t850_train.isel(time=0).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def create_training_data(data, lead_time_h, return_valid_time=False):\n",
    "    \"\"\"Function to split input and output by lead time.\"\"\"\n",
    "    X = data.isel(time=slice(0, -lead_time_h))\n",
    "    y = data.isel(time=slice(lead_time_h, None))\n",
    "    valid_time = y.time\n",
    "    if return_valid_time:\n",
    "        return X.values.reshape(-1, nlat*nlon), y.values.reshape(-1, nlat*nlon), valid_time\n",
    "    else:\n",
    "        return X.values.reshape(-1, nlat*nlon), y.values.reshape(-1, nlat*nlon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Train linear regression\n",
    "\n",
    "Now let's train the model. We will use scikit-learn for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def train_lr(lead_time_h, vars=['z', 't']):\n",
    "    \"\"\"Create data, train a linear regression and return the predictions.\"\"\"\n",
    "    X_train, y_train, X_test, y_test = [], [], [], []\n",
    "    if 'z' in vars:\n",
    "        X, y = create_training_data(data_z500_train, lead_time_h)\n",
    "        X_train.append(X); y_train.append(y)\n",
    "        X, y, valid_time = create_training_data(data_z500_test, lead_time_h, return_valid_time=True)\n",
    "        X_test.append(X); y_test.append(y)\n",
    "    if 't' in vars:\n",
    "        X, y = create_training_data(data_t850_train, lead_time_h)\n",
    "        X_train.append(X); y_train.append(y)\n",
    "        X, y, valid_time = create_training_data(data_t850_test, lead_time_h, return_valid_time=True)\n",
    "        X_test.append(X); y_test.append(y)\n",
    "    X_train, y_train, X_test, y_test = [np.concatenate(d, 1) for d in [X_train, y_train, X_test, y_test]]\n",
    "        \n",
    "    lr = LinearRegression()\n",
    "    lr.fit(X_train, y_train)\n",
    "    \n",
    "    mse_train = mean_squared_error(y_train, lr.predict(X_train))\n",
    "    mse_test = mean_squared_error(y_test, lr.predict(X_test))\n",
    "    print(f'Train MSE = {mse_train}'); print(f'Test MSE = {mse_test}')\n",
    "    preds = lr.predict(X_test).reshape((-1, len(vars), nlat, nlon))\n",
    "    \n",
    "    fcs = []\n",
    "    for i, v in enumerate(vars):\n",
    "        if v == 'z': data_mean = z500_mean; data_std = z500_std\n",
    "        else: data_mean = t850_mean; data_std = t850_std\n",
    "        fc = xr.DataArray(\n",
    "            preds[:, i] * data_std + data_mean, \n",
    "            dims=['time', 'lat', 'lon'],\n",
    "            coords={\n",
    "                'time': valid_time,\n",
    "                'lat': z500.lat,\n",
    "                'lon': z500.lon\n",
    "            },\n",
    "            name=v\n",
    "        )\n",
    "        fcs.append(fc)\n",
    "    return fcs, lr   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 3 days\n",
    "\n",
    "Here we train a model to directly predict the fields at 3 days lead time. Let's train a model that only predicts z or t and then a combined model. As we can see below, the model trained only on Z500 performs better than the combined model. But the same is not the case for T850. For the paper, we will use the combined model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train MSE = 0.43062442541122437\n",
      "Test MSE = 0.5234634280204773\n"
     ]
    }
   ],
   "source": [
    "# Train only with z\n",
    "fc_3d, lr_3d = train_lr(3*24, vars=['z'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(692.76963844)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_weighted_rmse(fc_3d[0], z500_valid).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train MSE = 0.3241349160671234\n",
      "Test MSE = 0.40045636892318726\n"
     ]
    }
   ],
   "source": [
    "# Train only with t\n",
    "fc_3d, lr_3d = train_lr(3*24, vars=['t'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(3.18926955)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_weighted_rmse(fc_3d[0], t850_valid).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train MSE = 0.31796592473983765\n",
      "Test MSE = 0.48072150349617004\n"
     ]
    }
   ],
   "source": [
    "# Train with z and t\n",
    "(fc_z500_3d, fc_t850_3d), lr_3d = train_lr(3*24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(713.78158629)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_weighted_rmse(fc_z500_3d, z500_valid).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(3.18461606)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_weighted_rmse(fc_t850_3d, t850_valid).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "xr.Dataset({'z': fc_z500_3d, 't': fc_t850_3d}).to_netcdf(f'{PREDDIR}/fc_lr_3d.nc');\n",
    "to_pickle(lr_3d, f'{PREDDIR}/saved_models/lr_3d.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### 5 days\n",
    "\n",
    "Combined model for 5 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train MSE = 0.3933986723423004\n",
      "Test MSE = 0.6094844341278076\n"
     ]
    }
   ],
   "source": [
    "# Train with z and t\n",
    "(fc_z500_5d, fc_t850_5d), lr_5d = train_lr(5*24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(813.53949079)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_weighted_rmse(fc_z500_5d, z500_valid).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(3.51787942)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_weighted_rmse(fc_t850_5d, t850_valid).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "xr.Dataset({'z': fc_z500_5d, 't': fc_t850_5d}).to_netcdf(f'{PREDDIR}/fc_lr_5d.nc');\n",
    "to_pickle(lr_3d, f'{PREDDIR}/saved_models/lr_5d.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Iterative forecast\n",
    "\n",
    "Finally, an iterative forecast. First, we train a model for 6 hours lead time and then construct an iterative forecast up to 120 hours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train MSE = 0.027540646493434906\n",
      "Test MSE = 0.03457337245345116\n"
     ]
    }
   ],
   "source": [
    "_, lr_6h = train_lr(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "to_pickle(lr_6h, f'{PREDDIR}/saved_models/lr_6h.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'z' (time: 17520, lat: 32, lon: 64)>\n",
       "dask.array<truediv, shape=(17520, 32, 64), dtype=float32, chunksize=(8760, 32, 64)>\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -87.19 -81.56 -75.94 -70.31 ... 75.94 81.56 87.19\n",
       "  * lon      (lon) float64 0.0 5.625 11.25 16.88 ... 337.5 343.1 348.8 354.4\n",
       "  * time     (time) datetime64[ns] 2017-01-01 ... 2018-12-31T23:00:00"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_z500_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "state = np.concatenate([data_z500_test.values.reshape(-1, nlat*nlon), \n",
    "                        data_t850_test.values.reshape(-1, nlat*nlon)], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(17520, 4096)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def create_iterative_fc(state, model, lead_time_h=6, max_lead_time_h=5*24):\n",
    "    max_fc_steps = max_lead_time_h // lead_time_h\n",
    "    fcs_z500, fcs_t850 = [], []\n",
    "    for fc_step in range(max_fc_steps):\n",
    "        state = model.predict(state)\n",
    "        fc_z500 = state[:, :nlat*nlon].copy() * z500_std + z500_mean\n",
    "        fc_t850 = state[:, nlat*nlon:].copy() * t850_std + t850_mean\n",
    "        fc_z500 = fc_z500.reshape((-1, nlat, nlon))\n",
    "        fc_t850 = fc_t850.reshape((-1, nlat, nlon))\n",
    "        fcs_z500.append(fc_z500); fcs_t850.append(fc_t850)\n",
    "\n",
    "    return [xr.DataArray(\n",
    "        np.array(fcs), \n",
    "        dims=['lead_time', 'time', 'lat', 'lon'],\n",
    "        coords={\n",
    "            'lead_time': np.arange(lead_time_h, max_lead_time_h + lead_time_h, lead_time_h),\n",
    "            'time': data_z500_test.time,\n",
    "            'lat': data_z500_test.lat,\n",
    "            'lon': data_z500_test.lon\n",
    "        }\n",
    "    ) for fcs in [fcs_z500, fcs_t850]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "fc_z500_6h_iter, fc_t850_6h_iter = create_iterative_fc(state, lr_6h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "xr.Dataset({'z': fc_z500_6h_iter, 't': fc_t850_6h_iter}).to_netcdf(f'{PREDDIR}/fc_lr_6h_iter.nc');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "rmses_z500_6h_iter = evaluate_iterative_forecast(fc_z500_6h_iter, z500_valid)\n",
    "rmses_t850_6h_iter = evaluate_iterative_forecast(fc_t850_6h_iter, t850_valid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEMCAYAAADXiYGSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXiU5b3/8XcyIXtC9mRCAiEIOIACBqVKBcWFHg0EWzU1rq3SugHq8aqRtuCRosZ6KrSIUrX+Dq3HtlYBjda4gBVcWBSVGGRJAgQy2SYJ2beZ5/cHNgdaEgfIrPm8rsvLJPdMnu+Xycxn5n6WO8AwDAMREZFvEejpAkRExDcoMERExCkKDBERcYoCQ0REnKLAEBERpygwRETEKQoMERFxSpCnC3C1hoZWHA7/ONUkPj4Sm63F02UMOH/sSz35Dn/s63R6CgwMIDY24oRjfh8YDofhN4EB+FUvx/LHvtST7/DHvlzRk6akRETEKQoMERFxigJDREScosAQERGnKDBERMQpCgwREXGK3x9WKyLib1o6e7A2dWBt6sR65Jv/N3X0/iw1Joz/yZs04Nt1W2Bs3LiRFStWYBgGDoeD+fPnc/nll1NeXk5+fj6NjY3ExMRQUFBARkYGQL9jIiL+yDAMmjqOCYTjguHo182dPcfdJyQoEHN0CCnRoViSo5gxLtkltQW4Y8U9wzA477zzePHFFxkzZgxff/011113HZ9++im33HILP/jBD8jJyWH9+vW88sorrFmzBoCbbrqpzzFn2WwtfnNSTmJiFLW1zZ4uY8D5Y1/qyXd4oi/DMGho76aioZ2KxnYqGjs41Pt1Oy2d9uNuHz7EREp0CKlDQ0mJOvp/c3Qo5ugQzENDiQ0bQkBAwID0FBgYQHx85AnH3PYJIzAwkObmow00NzeTlJREQ0MDJSUlvPDCCwBkZ2ezdOlS6uvrMQyjz7G4uDh3lS0ickoMw8DW2kVFYwcVje0camz/JiA6ONTYTmvX/4VCYACYo0NJjwlj/JlRpMWEYR4aSuo3nxqGhgYdFwie4pbACAgIYPny5dx5552Eh4fT2trK6tWrsVqtJCcnYzKZADCZTCQlJWG1WjEMo8+xkwmMvpLSVyUmRnm6BJfwx77Uk+84nb46e+zsr2tjX00Le2ua2VvTQlltKwdsrbQdEwqmwADSY8MYER/B+WckMCI+nIz4CEbEh5MWG05w0MAeg+SKx8otgdHT08Pq1atZtWoVWVlZfPrpp9x77708/vjjLt+2pqS8nz/2pZ58h7N9dXTb2V/fRnl9G+W2//vvUGM79m9eYgKA1KGhZMSFM3FCCukxoaTFhDE8NoyUqBCCTCcKBYMjDa0e6elEPD4ltWvXLmpqasjKygIgKyuLsLAwQkJCqK6uxm63YzKZsNvt1NTUYDabMQyjzzEREVdp67JTWtd6fDDUt2E90sE/33qaAiA9NozMhAguGZPAyPgIRsaHMyI2jNAhJo/W70puCYyUlBSqqqooKysjMzOT0tJS6urqGDFiBBaLhcLCQnJycigsLMRisfROOfU3JiJyuupauyiur2Hb3lr21LSyp7aFiob23mAINgUwIi6cCSlRZI9PJjM+nJHx4aTHhDHkhJ8W/JtbjpICeO2113j22Wd7d9wsWLCASy+9lNLSUvLz82lqaiI6OpqCggIyMzMB+h1zlqakvJ8/9qWevIvdYVDR2M6emhb21Layu6aFPTUt1Ld1994mdWgoY5MiGZMYwejECDLjI0gdGoop0PM7m0+Wq6ak3BYYnqLA8H7+2Jd68pzOHgf76v4vFPbUtLKvroX2bgcAQYEBZMaHMyYpkrFJkZw3OpHEIYFEhfrPecw+vQ9DRMQVDMPg8JEOiq3NFFubKLY2s7umhZ5v3iRGBJsYkxRJzllmxiRGMCYpksz48OOmk3wlCL2BAkNEfEZLZw9fWZsprmr6JiSaaWw/Oq0UNiQQS3IUeVlpjE+JZGxyJKnRoV5x/oK/UGCIiFfqcRiU1bVSXNVMceXRgNhf39a7Q3pkXDgXZsYxITWaCSlRZCZEEOSD+xt8iQJDRLxCR7edLyub+PTQET4/dIRd1c29+x1iwoYwwRzFLEsiE1KiGZcS5Vf7HHyF/sVFxCM6uu3stDbxacURPq1opNjaTI/DwBQAY5IimTMhhQnmaCaYoxg2VFNL3kCBISJu0dnjYGdlE59WNPLpoSMUW5vothsEBsCZyVHkZQ3jnPQYJg2LJiJYL03eSI+KiLhEZ4+DYus3AVFxNCC6vgmIsUmR5E4expT0GCYOiyYyRC9FvkCPkogMCMMw2Fvbyofl9Ww50MDOyqMBEQCcmRzJNZOGkZU+lMlpQxUQPkqPmoicsrYuO1sPNPBheT0flddT09IFwJjECK6elEpWegyThw3VDmo/oUdRRE7KwYZ2NpfZ+Ki8ns8OHaHbbhARbGLqiFimjYzjgpGxJESGeLpMcQEFhoj0q6vHwWeHGvnsk4O8+1UVFY0dAGTEhXHtpGFMy4xl0rChg/JifIONAkNE/k11c+fRaaayerYebKC920FwUCBZaUP54TnDuGBkHGkxYZ4uU9xMgSEiAFQ1dfDO7lre2V3LruoWAJKjQrhiXDLTRsbxH+ek03KkzcNViicpMEQGsbqWTt7dU8c7u2v5srIJAEtyJHd9N4PvjopnVHx47wlzYcEmWjxZrHicAkNkkGlo62LD3qMh8VnFEQzgjIQI7piWwWVjE0mP1VSTnJgCQ2QQaOro5v29Nt7ZXcu2gw3YDRgRG8at3xnOZWcmkhkf4ekSxQcoMET8VEtnDx+UHg2JT/Y30OMwSB0ayg3npnP52ERGJ0bo+kxyUhQYIn6kx+5gc1k9b+6q4cMyG112g6TIYHInD+OyMxMZlxypkJBTpsAQ8QOHj7Tz2s4qXiuupq61i7jwIcw9y8zlZyZyVmo0gQoJGQAKDBEf1WN38EGpjbU7q9iyv4GAALhgZBxzzzIzLTNOiwnJgHNLYBw6dIi77rqr9/vm5mZaWlrYunUr5eXl5Ofn09jYSExMDAUFBWRkZAD0OyYyWB1qbGfdzipeL66ivq2bpMhgbjt/OHMmpJASHerp8sSPuSUw0tLSWL9+fe/3y5Ytw263A7BkyRLy8vLIyclh/fr1LF68mDVr1nzrmMhg0m138I99NtbttLLlQCOBATBtZBxXnW3m/JH6NCHu4fYpqa6uLl5//XWef/55bDYbJSUlvPDCCwBkZ2ezdOlS6uvrMQyjz7G4uDh3ly3iERUN7azbaeX14moa2rtJiQrhpxeMYPaEFJKjdIE/cS+3B8aGDRtITk5m/PjxFBcXk5ycjMlkAsBkMpGUlITVasUwjD7HFBjiz7rtDjburWPtziq2H2zEFAAXjopn7tlmvjMiFpM+TYiHuD0wXnnlFX7wgx+4bXvx8ZFu25Y7JCZGeboEl/DHvk62p+aObv53y0H+8GE51U2dDIsJ4/7Lx3DNlHSSvWTfhD8+TuCffbmiJ7cGRnV1Ndu2bePxxx8HwGw2U11djd1ux2QyYbfbqampwWw2YxhGn2Mnw2ZrweEwXNGO2yUmRlFb2+zpMgacP/Z1Mj3VtXTy5x2VvPJFJS2ddqYMj+HBS0dzfkbs0cNhO7upre12ccXfzh8fJ/DPvk6np8DAgD7faLs1MNauXcuMGTOIjY0FID4+HovFQmFhITk5ORQWFmKxWHqnnPobE/F1Bxva+eO2Ct4oqcbuMJg5OoEbz01nXIr/vdsV/+D2wPj5z39+3M8eeugh8vPzWbVqFdHR0RQUFDg1JuKrvqpqZs3WCjburWOIKYDZ41O4YUqaLvonXi/AMAz/mK/pg6akvJ8/9vWvPRmGwScHGliztYLtFUeIDDFx9cRUcs8ZRkJEsAcrdZ4/Pk7gn335xZSUyGDT4zB4b3cta7ZVsKe2lcTIYBZMH8lVZ5uJDNHTT3yL/mJFXKC9y85fdxzmxe2HqGzqJCMujF9ePobvWZIIDtLa1+KbFBgiA6ij286fPzvMSzsqqW/t4ixzNPddPIoLR8XrAoDi8xQYIgPA7jAo/KqK3390gJqWLi4am0jepFQmDYvW5cTFbygwRE6DYRhsLqvnd5vKKbe1McEcxdIrz2TW5HS/25EqosAQOUXF1iZ++0E5Ow4dYXhsGI/NtjBzdII+UYjfUmCInKSDDe2s2lzOe3vqiAsfws8uOYOrzkohyKSd2eLfFBgiTrK1dvHcxwdYu7OKYFMA884fzvVT0ogI1tNIBgf9pYt8i7YuOy9uP8Sfth+is8fO3LPN3Hb+CJ854U5koCgwRPrQY3ewbmcVz358gPq2bmaOTuDO72YwIi7c06WJeIQCQ+RfGIbBxn02ntpUzsGGdiYPi+aJnPGclRrt6dJEPEqBIXKM/bY2Hnl3LzsOHWFkfDj/PXc8F2bG6cgnERQYIgB09Tj4n60VvLD1IGFDTDx42WjmTEjRWtkix1BgyKC349ARHnlnD/vr25l1ZiL3XTyKuHDt0Bb5VwoMGbSaO3r43aYy1n5ZRWp0CCu+P4ELRmqBLpG+KDBk0DEMg/f21PHExlIa2rq4PiuNn04bQdgQk6dLE/FqCgwZVKqaOih4bx+by+o5MymS5VeN58xkLYkq4gwFhgwKdofBXz+v5OnN5RgG3DMjk9xzhmmntshJUGCI39td08Ij7+ylpKqZC0bG8sAlo0kdGurpskR8jgJD/FZHt51nPz7Ai9sPMTRsCMuuPJPLxibqnAqRU6TAEL/0yf56Hn13H5VHOsg5K4X5F45kaNgQT5cl4tPcFhidnZ088sgjfPzxx4SEhDBp0iSWLl1KeXk5+fn5NDY2EhMTQ0FBARkZGQD9jomcSFuXnSc27OP1r6oZHhvGM9eeTVZ6jKfLEvELbguMX//614SEhFBUVERAQAB1dXUALFmyhLy8PHJycli/fj2LFy9mzZo13zom8q/21LSwqHAXBxva+dHUdG79zghCgrRGhchAccuzqbW1lXXr1rFw4cLe+eOEhARsNhslJSVkZ2cDkJ2dTUlJCfX19f2OiRzLMAxe/rySH/3vDlq77Ky65mzu/O5IhYXIAHPLJ4yKigpiYmJYuXIlW7ZsISIigoULFxIaGkpycjIm09ETpkwmE0lJSVitVgzD6HMsLs75s3Hj4yNd0pOnJCb65zkDp9rXkbZu8l/9kr8XV3HR2ESeuGYiCZEhA1zdqfHHx8ofewL/7MsVPbklMHp6eqioqGDcuHE88MADfPHFF9x+++2sWLHC5du22VpwOAyXb8cdEhOjqK1t9nQZA+5U+/qysolfvLGLmpYuFkwfyfVT0jDau6ht73JBlSfHHx8rf+wJ/LOv0+kpMDCgzzfabgmM1NRUgoKCeqeXJk6cSGxsLKGhoVRXV2O32zGZTNjtdmpqajCbzRiG0eeYDG4Ow+CP2w7x9OZykqNCeO6HE5lg1loVIq7mlkneuLg4pk6dyocffggcPfrJZrORkZGBxWKhsLAQgMLCQiwWC3FxccTHx/c5JoOXrbWLha8Us3JTORePTuBPN2YpLETcJMAwDLfM11RUVLBo0SIaGxsJCgrinnvuYcaMGZSWlpKfn09TUxPR0dEUFBSQmZkJ0O+YszQl5f2c7WvrgQYW/303LZ093HdRJledbfbak/D88bHyx57AP/ty1ZSU2wLDUxQY3u/b+upxGDz70X5e2FJBRlw4j2RbOCMxwo0Vnjx/fKz8sSfwz758eh+GyKmqaurgF298zReVTcyZkMz9M8/QZchFPESBIV7rH/vqeLhoDz12g6VXnMn3LEmeLklkUFNgiNfp6nHw2w/K+MuOSs5MimRZtoXhsWGeLktk0FNgiFepa+3i/nVf8VVVMz88ZxjzLxxJsM7YFvEKCgzxGqV1rdy7tpiGtm4K5oxj5ugET5ckIsdQYIhX2HKggQdeKyF0iInVuRMZl+J/l2oQ8XUKDPG4P289yM/XFTMyLpwnrxpPSrRWwxPxRgoM8RiHYfDUpv2s2VbBdzJieTTbQmSI/iRFvJWeneIRHd12HnprN+/tqSNv6nDmXzCCoEDvPGtbRI5SYIjb1bd18Z/rvuIrazMLZ2Ryz/fOpK6uxdNlici3UGCIW5XZWrn31WJs3xwJdfHoBK+9HpSIHE+BIW6z9UADD7xeQrApkNW5ExmvI6FEfIoCQ9zitZ1VPPLuXkbEhrH8+xMw60goEZ+jwBCXchgGT2/ez//bWsF3RsTy6GwdCSXiq/TMFZfp6LbzX2/t4d09tVx1dgo/m3kGQSZd5kPEVykwxCXq245eE2qntZkF00dyw5Q07dwW8XEKDBlw+21tLFxbjK21i4LZFmaOSfR0SSIyABQYMqD21bZyx8tfEhgAq689m/Fab1vEbygwZMD8MyyCTQE8fe1ErWEh4mec3gNpGAZ//etfuemmm5g9ezYA27Zt480333RZceI7FBYi/s/pwFixYgV/+9vfyM3NxWq1ApCSksJzzz3n1P1nzpzJ9773PXJycsjJyWHTpk0AlJeXk5uby6xZs8jNzWX//v299+lvTLyHwkJkcHA6MNauXcszzzzDlVde2Xu0S1paGhUVFU5v7Le//S3r169n/fr1XHjhhQAsWbKEvLw8ioqKyMvLY/Hixb23729MvIPCQmTwcDow7HY7ERERAL2B0draSnh4+Clv3GazUVJSQnZ2NgDZ2dmUlJRQX1/f75h4B4WFyODidGDMmDGDRx99lK6uLuDoPo0VK1Zw8cUXO72x+++/n9mzZ/PQQw/R1NSE1WolOTkZk8kEgMlkIikpCavV2u+YeJ7CQmTwcfooqQcffJCf/exnZGVl0dPTw+TJk5k2bRoFBQVO3f/FF1/EbDbT1dXFsmXLePjhh7nllltOtW6nxcdHunwb7pSY6PkL9u2yNnHXKzsJHWLipZ98h5EJEaf9O72hr4GmnnyHP/blip6cDozIyEhWrVqFzWbj8OHDmM1mEhOdPyHLbDYDEBwcTF5eHnfccQcPPvgg1dXV2O12TCYTdrudmpoazGYzhmH0OXYybLYWHA7jpO7jrRITo6itbfZoDXtrW7jz5Z0EmwJYdfVZRBqO067JG/oaaOrJd/hjX6fTU2BgQJ9vtJ2ekqqvr6e1tZX4+HjGjx/PBx98wLp163A4HN9637a2NpqbjxZvGAZvvvkmFouF+Ph4LBYLhYWFABQWFmKxWIiLi+t3TDzj2LB45tqJpGsaSmRQCTAMw6m339dccw3/9V//xbhx43jiiSfYuHEjQUFBTJ06lUWLFvV734qKCubPn4/dbsfhcDBq1Ch+8YtfkJSURGlpKfn5+TQ1NREdHU1BQQGZmZkA/Y45S58wBoYrw0Lv8HyDP/YE/tmXqz5hOB0Y5557Llu3biUgIIDp06fz5z//mfDwcLKzs9m8efMpFeYOCozT5+pPFnrC+gZ/7An8sy9XBYbT+zACAwPp7u6mvLycqKgoUlNTcTgctLa2nlJR4hs0DSUi/+R0YEyfPp2FCxfS2NjIFVdcAcC+fftITk52WXHiWQoLETmW04GxbNky1q5dS1BQEDk5OQA0NDQwf/58lxUnnqOwEJF/5XRgBAcHk5ube9zPpk6dOuAFiecpLETkRJwOjObmZtasWcOuXbtoa2s7buwPf/jDgBcmnlHV1MH8V4oVFiLyb5wOjIULF2K327nssssICQlxZU3iIS2dPdyztpjOHjvPXzdJYSEix3E6MD7//HO2bNnCkCFDXFmPeEiPw2BR4S7217ez4vsTyIw//ct9iIh/cfpM76ysLEpLS11Zi3iIYRg8sWEfH+9vIP+SM5g6ItbTJYmIF3L6E8Zjjz3GvHnzmDhxIvHx8ceN3X333QNemLjPS58d5pUvrNx0bhpzzz65a3WJyODhdGA8+eSTVFVVkZaWRktLS+/P/7k2hvimf+yrY/n7ZVw8OoG7Lhzp6XJExIs5HRhvvPEGRUVFJCUlubIecaNd1c384o2vsaRE8fB/jCVQ4S8i/XB6H0Z6ejpBQU7ni3i56uZO7lv7FTFhQ/jvueMJHWLydEki4uWcToCcnBzuvPNObrjhhn/bh3H++ecPeGHiOq1dPdy7tpj2bjvPXTeJhIhgT5ckIj7A6cBYs2YNJpOJ3/zmN8f9PCAggPfee2/ACxPX6HEY/OKNrymra+XJ70/gjAFYLU9EBgenAsNut1NfX8+nn35KcLDejfqy5e+XsrmsnvxLz+D8DC1GJSLOc2ofhslkYuTIkTQ0NLi6HnGhv3x2mL/sqCQvaxg/mJjq6XJExMc4PSU1e/Zsbr/9dm666SZSUlKOG9M+DO+3uczGb94vZcaoeBZMP7lVC0VE4CQC46WXXgLgd7/73XE/1z4M77e7poVFhbsYmxTJ0ivPxBSow2dF5OQ5HRgbNmxwZR3iIjXNndy3tpiokCD+e+54wnT4rIicIp1Y4cfauuzct+4rWjrtPHfdRBIjdZVhETl1Cgw/ZXcY/OKNXeytbeE3cycwOvHEi7qLiDjL6TO9B8rKlSsZO3Yse/bsAaC8vJzc3FxmzZpFbm4u+/fv771tf2PSvxX/KGNTWT3/efEZTMvU4bMicvrcGhhfffUVn3/+Oamp/3dI55IlS8jLy6OoqIi8vDwWL17s1Jj07W+fV/LSZ4fJnZzKtZN1+KyIDAy3BUZXVxcPP/wwS5Ys6b3Crc1mo6SkhOzsbACys7MpKSmhvr6+3zHp21dVzfz3xlKmjYzj3otGebocEfEjbtuHsWLFCubMmUN6enrvz6xWK8nJyZhMR4/cMZlMJCUlYbVaMQyjz7G4OOenWOLj/WvuPjExqs+xI+3d/PLv20iODuWpG7OICfeds/L768tXqSff4Y99uaIntwTGjh072LlzJ/fff787Nnccm60Fh8Nw+3ZdITExitra5hOOGYbBg4W7qGzs4Pe5E+lu7aS2tdPNFZ6a/vryVerJd/hjX6fTU2BgQJ9vtN0yJbVt2zbKysq45JJLmDlzJlVVVdx6660cPHiQ6upq7HY7cPSaVTU1NZjNZsxmc59j8u/+9oWV9/bUcdd3Mzg7NdrT5YiIH3JLYPzkJz9h8+bNbNiwgQ0bNpCSksLzzz/PFVdcgcViobCwEIDCwkIsFgtxcXHEx8f3OSbH213TwpPvl3LByFiun5Lm6XJExE95/DyMhx56iPz8fFatWkV0dDQFBQVOjclRrV09LCrcRUzYEB76nlbNExHX8UhgHHuZkVGjRvHyyy+f8Hb9jcnR/RaPvrOXQ43tPH3t2cT60E5uEfE9bj9xTwbOa8VVFH1dy08uGME5aTGeLkdE/JwCw0ftq2vl1xtKOXd4DLecN9zT5YjIIKDA8EHt3XYWvb6LiGATD1+hy5WLiHt4fKe3nLxfv7eP/fVt/O7qs0iI0H4LEXEPfcLwMW+WVPP6V9X86DvDmToi1tPliMggosDwIftqWnjs3b1MThvKvPNHeLocERlkNCXlIzq67dz91y8JNgXyqyvOJEj7LUTEzRQYPuLJ98v4uqqZ5d+fQFKUVs4TEffTlJQPePvrGl790spPZ2QybaQujSIinqHA8HKHGtt55J29nGWO4v7Lx3q6HBEZxBQYXqyrx8GDr+8iMCCAZdkWhpj0cImI5+gVyIv99oMyvq5pYfGsMZijQz1djogMcgoML/X+3jr+sqOSH54zjItGJ3i6HBERBYY3qjzSwcNFe7AkRzL/wpGeLkdEBFBgeB2HYbD4za9xGAaPZFsIDtJDJCLeQa9GXmbdziq+qGziPy8eRVpMmKfLERHppcDwInWtXaz8oJys9KFkj0/2dDkiIsdRYHiR5e+X0tFjJ//S0QRoqVUR8TIKDC/x8f56ir6u5Zbz0smIC/d0OSIi/0aB4QU6uu0UvLuP4bFh3KzV80TES7nt4oN33nknhw4dIjAwkPDwcH75y19isVgoLy8nPz+fxsZGYmJiKCgoICMjA6DfMX/yhy0HOXykg6evOZsQHRUlIl7Kba9OBQUFvPbaa6xbt44f//jHLFq0CIAlS5aQl5dHUVEReXl5LF68uPc+/Y35i9K6VtZsO8SV45OZMjzG0+WIiPTJbYERFRXV+3VLSwsBAQHYbDZKSkrIzs4GIDs7m5KSEurr6/sd8xcOw+DRd/YSGWzinumZni5HRKRfbl0P4+c//zkffvghhmHw3HPPYbVaSU5OxmQyAWAymUhKSsJqtWIYRp9jcXH+cYnv9d+cc7F41hhiwod4uhwRkX65NTCWLVsGwLp163j88cdZuHChy7cZHx/p8m2citrmTlZuKmfqyDh+dNEZTh9Gm5gY9e038kH+2Jd68h3+2JcrevLIintz585l8eLFpKSkUF1djd1ux2QyYbfbqampwWw2YxhGn2Mnw2ZrweEwXNTJqfvlm1/T3m3n/hmZ1NW1OHWfxMQoamubXVyZ+/ljX+rJd/hjX6fTU2BgQJ9vtN2yD6O1tRWr1dr7/YYNGxg6dCjx8fFYLBYKCwsBKCwsxGKxEBcX1++Yr9uyv4G3dtVw87npZMTrnAsR8Q1u+YTR3t7OwoULaW9vJzAwkKFDh/LMM88QEBDAQw89RH5+PqtWrSI6OpqCgoLe+/U35qs6uu089t5ehseGcctUnXMhIr4jwDAM75uvGUDeNiX19OZy/rClglXXnMW5w2NP6r7++NEZ/LMv9eQ7/LEvn56SkqPKbN+cczEu6aTDQkTE0xQYbvLPcy4igk0snKFzLkTE9ygw3OT14io+P9zEghmZxIYHe7ocEZGTpsBwg/q2Ln77QTmT04YyW+tciIiPUmC4wfL3y2jrsrNI61yIiA9TYLjYlgMN/H1XDTefp3MuRMS3KTBcqLPHQcG7e0mPCeVHOudCRHycAsOFXthykIrGDvIvHa11LkTE5+lVzEXKbW38z9YKrhiXxHkjdM6FiPg+BYYLHD3nYg8RwSbu0TkXIuInFBguUFhczY7DTSyYrnMuRMR/KDAGWGtXDys3lTNpWDTZE3TOhYj4DwXGAHtx+yEa2ru556JRBOqcCxHxIwqMAVTf1sWL2w9zyZgExqf43wpeIjK4KTAG0B8+OUhnj53bp2V4uhQRkQGnwBggh4+089xKYTgAAAxpSURBVMoXVmZPSCEjTmd0i4j/UWAMkN9/dABTYADzzh/h6VJERFxCgTEA9ta28PeSGnInDyMpKsTT5YiIuIQCYwCs2ryfyJAgbj4vzdOliIi4jALjNH1+6Aiby+q5+bx0okOHeLocERGXUWCcBsMw+N2mchIjg8mdnOrpckREXMotgdHQ0MC8efOYNWsWs2fP5u6776a+vh6A8vJycnNzmTVrFrm5uezfv7/3fv2NeYNNZfV8WdnEbeePIHSIydPliIi4lFsCIyAggNtuu42ioiJef/110tPTeeKJJwBYsmQJeXl5FBUVkZeXx+LFi3vv19+Yp9kdBk9tKmd4bBhzJqR4uhwREZdzS2DExMQwderU3u8nTZpEZWUlNpuNkpISsrOzAcjOzqakpIT6+vp+x7zB33dVU2Zr487vZhAUqEuAiIj/C3L3Bh0OBy+99BIzZ87EarWSnJyMyXR0OsdkMpGUlITVasUwjD7H4uLinN5efHzkgPfQ2WPnuU8qODttKLkXjHTrOt2Jif55yRF/7Es9+Q5/7MsVPbk9MJYuXUp4eDg33HADJSUlLt+ezdaCw2EM6O/8308PcbixnUWXnkFdXcuA/u7+JCZGUVvb7LbtuYs/9qWefIc/9nU6PQUGBvT5RtutgVFQUMCBAwd45plnCAwMxGw2U11djd1ux2QyYbfbqampwWw2YxhGn2Oe1NLZwwtbKpg6IkYr6YnIoOK2w2qffPJJiouLeeqppwgOPrqoUHx8PBaLhcLCQgAKCwuxWCzExcX1O+ZJL24/RGN7N3ddONKjdYiIuFuAYRgDO19zAnv37iU7O5uMjAxCQ0MBSEtL46mnnqK0tJT8/HyampqIjo6moKCAzMyjy5r2N+asgZySsrV2cdXzW5k2Mp5HZ1sG5HeeDH/86Az+2Zd68h3+2JdPT0mNHj2a3bt3n3Bs1KhRvPzyyyc95gl/+OQgXT0O7vhuhqdLERFxO53p7aRDje28+qWVnLPMDI8N83Q5IiJup8Bw0upvLl9+2/nDPV2KiIhHKDCcsKemhaJdNfzwnGEkRury5SIyOCkwnLBq836iQoO4+dx0T5ciIuIxCoxv8dmhRj4sr+eW89KJCnX7eY4iIl5DgdEPwzBY+UE5SZHBXDNJly8XkcFNgdGPD0pt7LQ2M0+XLxcRUWD0xe4weGrzfkbEhpGty5eLiCgw+vJGSTXluny5iEgvBcYJdPY4+P1HBxifEsXFoxM8XY6IiFdQYJzA3z6vpLq5k7svdO9aFyIi3kyBcQKfVjRyYWYcU4bHeLoUERGvoRMLTuCRbAtBJmWpiMixFBgnoENoRUT+nd5Gi4iIUxQYIiLiFAWGiIg4RYEhIiJOUWCIiIhTFBgiIuIUvz+sNtDPrgPlb/38kz/2pZ58hz/2dao99Xe/AMMwjFMtSEREBg9NSYmIiFMUGCIi4hQFhoiIOEWBISIiTlFgiIiIUxQYIiLiFAWGiIg4RYEhIiJOUWCIiIhTFBheqKGhgXnz5jFr1ixmz57N3XffTX19PQDl5eXk5uYya9YscnNz2b9/v2eLPQUrV65k7Nix7NmzB/D9njo7O1myZAmXX345s2fP5pe//CXg231t3LiRuXPnkpOTw+zZs3n77bcB3+upoKCAmTNnHvf3Bv334e09nqin/l4zYAB7MsTrNDQ0GJ988knv94899pjx4IMPGoZhGDfeeKOxbt06wzAMY926dcaNN97okRpPVXFxsXHrrbcaF110kbF7927DMHy/p6VLlxrLli0zHA6HYRiGUVtbaxiG7/blcDiMKVOm9D4+u3btMiZNmmTY7Xaf62nbtm1GZWWlcfHFF/f2Yxj9Pzbe3uOJeurvNcMwBq4nBYYPeOutt4ybb77ZqKurM7Kysoyenh7DMAyjp6fHyMrKMmw2m4crdE5nZ6dx7bXXGgcPHuz9Y/f1nlpaWoysrCyjpaXluJ/7cl8Oh8M477zzjO3btxuGYRhbt241Lr/8cp/u6dgX1/768KUe/zUEj/XP1wzDGNi/Rb+/Wq2vczgcvPTSS8ycOROr1UpycjImkwkAk8lEUlISVquVuLg4D1f67VasWMGcOXNIT0/v/Zmv91RRUUFMTAwrV65ky5YtREREsHDhQkJDQ322r4CAAJYvX86dd95JeHg4ra2trF692ucfq3/qrw/DMHy+x2NfM2Bgn2Pah+Hlli5dSnh4ODfccIOnSzktO3bsYOfOneTl5Xm6lAHV09NDRUUF48aN49VXX+X+++9n/vz5tLW1ebq0U9bT08Pq1atZtWoVGzdu5Omnn+bee+/16Z4GE1e+ZigwvFhBQQEHDhxg+fLlBAYGYjabqa6uxm63A2C326mpqcFsNnu40m+3bds2ysrKuOSSS5g5cyZVVVXceuutHDx40Gd7AkhNTSUoKIjs7GwAJk6cSGxsLKGhoT7b165du6ipqSErKwuArKwswsLCCAkJ8dmejtXf88iXn2Pw768Z0H+/J0uB4aWefPJJiouLeeqppwgODgYgPj4ei8VCYWEhAIWFhVgsFp/4qPyTn/yEzZs3s2HDBjZs2EBKSgrPP/88V1xxhc/2BBAXF8fUqVP58MMPgaNHo9hsNjIyMny2r5SUFKqqqigrKwOgtLSUuro6RowY4bM9Hau/55EvP8dO9JoBA/u6oQWUvNDevXvJzs4mIyOD0NBQANLS0njqqacoLS0lPz+fpqYmoqOjKSgoIDMz08MVn7yZM2fyzDPPMGbMGJ/vqaKigkWLFtHY2EhQUBD33HMPM2bM8Om+XnvtNZ599lkCAo6uvrZgwQIuvfRSn+vpV7/6FW+//TZ1dXXExsYSExPDG2+80W8f3t7jiXpavnx5n68ZMHA9KTBERMQpmpISERGnKDBERMQpCgwREXGKAkNERJyiwBAREacoMERExCkKDBGOnhfy0UcfDejvfPXVV7nuuutO+n7bt29n1qxZA1qLyEBQYIh42NixYzlw4EDv91OmTKGoqMiDFYmcmAJDREScosAQOYbD4eD3v/89l156KVOnTmXhwoU0Njb2ji9YsIBp06aRlZXF9ddfz969e3vHGhoauP322znnnHO4+uqrOXjw4Ldu7/rrrwcgJyeHyZMn8+abb7JlyxamT5/ee5uZM2fy3HPPMXv2bCZNmsSiRYuoq6vjtttuY/Lkydxyyy0cOXKk9/aff/45P/zhD5kyZQpz5sxhy5YtA/FPI6LAEDnWmjVrePfdd/nTn/7Epk2bGDp0KA8//HDv+PTp0ykqKuLjjz9m3Lhx3H///b1jDz/8MCEhIWzevJlHHnmEV1555Vu39+KLLwKwfv16duzYwRVXXHHC27399tu88MILFBUVsXHjRubNm8d9993Hli1bcDgc/PGPfwSgurqan/70p9xxxx1s3bqVBx54gAULFhy3XKfIqVJgiBzjL3/5C/feey8pKSkEBwdz9913U1RURE9PDwBXX301kZGRBAcHM3/+fL7++muam5ux2+28/fbbLFiwgPDwcMaMGcNVV101YHXdcMMNJCQkkJyczJQpUzj77LMZN24cwcHBXHbZZZSUlABHg2f69OnMmDGDwMBApk2bxoQJE/jHP/4xYLXI4KUV90SOUVlZyV133dW7lgBAYGAgNpuNhIQEnnzySd566y3q6+t7b9PQ0EBHRwc9PT3HrTGQmpo6YHUlJCT0fh0SEnLc96Ghob2LG1VWVvLWW2+xcePG3vGenh6mTp06YLXI4KXAEDlGSkoKjzzySO/iQcdat24d7733Hi+88AJpaWk0Nzdz7rnnYhgGcXFxBAUFYbVaGTVqFHB0aUx3M5vN5OTk8Ktf/crt2xb/pykpkWNcd911LF++nMOHDwNQX1/Pu+++C0BrayvBwcHExsbS3t7Ob37zm977mUwmLrvsMlauXEl7ezv79u1j7dq1Tm0zISGBioqKAal/zpw5bNy4kU2bNmG32+ns7GTLli1UVVUNyO+XwU2BIXKMm266iZkzZ/LjH/+YyZMnc+211/Lll18CMHfuXFJTU7nwwgu58sormTRp0nH3Xbx4MW1tbUybNo38/Hy+//3vO7XNu+++m/z8fKZMmcKbb755WvWbzWZWrVrF6tWrOf/885kxYwbPP/88DofjtH6vCGgBJRERcZI+YYiIiFO001vExbZv3868efNOOLZjxw43VyNy6jQlJSIiTtGUlIiIOEWBISIiTlFgiIiIUxQYIiLiFAWGiIg45f8DcMomwvyU6jYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rmses_z500_6h_iter.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEcCAYAAAAydkhNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xUdb7/8VcmlRTSKwmdIEgoBkFERUIVgaArduxRr7R1l0XRq66oSLy7gFIuiuWn68r1chUWCFIirAIqCCJIDzWN9IT0NnN+f6BZsxAMkMxkJu/n45EHyZwzcz4fksw7p3zP18kwDAMREZELMNm6ABERabkUEiIi0iCFhIiINEghISIiDVJIiIhIgxQSIiLSIIWEOLy4uDi++eabZt3Gs88+y/z585t1GyK2oJAQaaEMw2D+/PnceOONxMbGMmnSJFJSUuqWT5o0iZiYGPr160e/fv0YNWpUved/++23jB49mj59+jBp0iQyMjKs3YI4AIWESAv1xRdf8Nlnn/HJJ5+wc+dO+vbty8yZM+ut8+KLL7Jnzx727NnDhg0b6h4vKChgypQpTJ8+nZ07d9KrVy+efvppa7cgDkAhIa2KxWLhnXfeYfjw4QwcOJDp06dTVFQEwKOPPsrHH39cb/3x48ezceNGAI4fP87DDz/MgAEDGDVqFOvWrWvWWtPT04mNjSUqKgpnZ2fGjx/PsWPHGvXcTZs20a1bN2655Rbc3d2ZOnUqhw8f5vjx481aszgehYS0Kh999BHJycl8/PHHbN26FV9fX2bPng3AuHHjWLt2bd26x44dIzMzk5tvvpny8nIeeeQRxo4dyzfffMO8efN4+eWX6x3+aciuXbvo379/gx+7du264PNuvfVWUlNTOXnyJDU1NaxcuZIbb7yx3jp//etfGThwIHfffTc7duyoezwlJYXu3bvXfe3p6Un79u0bHTIiv3CxdQEi1vTpp5/y4osvEhYWBsCUKVMYOnQotbW1DB8+nD//+c9kZGTQrl071qxZw4gRI3BzcyM5OZl27drxu9/9DoCrr76aUaNGsWHDBrp163bRbV4sCC4mODiY2NhYRo8ejbOzM2FhYXz44Yd1y2fMmEGXLl1wc3MjKSmJJ598kn/84x+0b9+e8vJyAgIC6r2et7c3ZWVll1yHtG4KCWlVMjMzmTx5MibTv3aiTSYT+fn5hIaGMmTIEJKSknj88cdJSkrilVdeASAjI4N9+/bRv3//uueZzWbGjx/fbLUuXryY/fv389VXXxEUFMTq1at58MEHSUpKok2bNvTp06du3dtuu421a9fy1VdfMWnSJDw9PSktLa33emVlZXh5eTVbveKYFBLSqoSFhTFnzhxiY2MvuHzs2LEsWrSIa6+9lsrKSgYOHAhAeHg41157LR988MElb3PXrl0kJCQ0uHzZsmX1wucXhw8f5pZbbqnb67n99tuZM2cOx44dIyYm5rz1nZyc+OWmzt26dWPlypV1y8rLy0lNTaVr166XXL+0bjonIa3KPffcw4IFC+ouBy0oKCA5Oblu+ZAhQ8jMzOStt95izJgxdXscN998M6dOnWLVqlXU1NRQU1PDvn37GnUiuH///nVXIF3o40IBARATE8P69evJy8vDYrGwatUqamtr6dChA8XFxWzdupWqqipqa2tZvXo1u3bt4oYbbgBgxIgRpKSksGHDBqqqqli8eDHdu3enS5cuV/pfKK2M9iSkVXnggQcwDINHHnmEnJwcAgMDGTNmDMOHDwfAzc2NESNG8Nlnn9W7ZNTb25v33nuPuXPnMnfuXAzDoHv37syaNavZak1ISCA/P58JEyZQXl5Ohw4deOutt2jbti0FBQUsWLCAEydO4OzsTOfOnVm8eDGdO3cGICAggIULFzJ79mz+9Kc/0adPH+bNm9dstYrjctKkQyIi0hAdbhIRkQYpJEREpEEKCRERaZBCQkREGqSQEBGRBikkRESkQQ43TqKwsAyLxXGu6g0M9CY/v/S3V7Qj6sl+OGJfjtgTXH5fJpMT/v4N367F4ULCYjEcKiQAh+sH1JM9ccS+HLEnaJ6+dLhJREQapJAQEZEGWe1w01NPPUV6ejomkwlPT09eeOEFevToUW+dhQsX8sknnxASEgLANddcw0svvWStEkVE5N9YLSQSExPx8fEBIDk5meeee67erYx/MWHCBJ555hlrlSUiIhdhtcNNvwQEQGlpKU5OTtbatIiIXCar3gX2+eefZ/v27RiGwbvvvnvetI8LFy5kxYoV+Pr6EhwczNSpU+nXr5+1yhMRsUuGYTTbH942uVX4qlWrSEpKYtmyZfUez83Nxc/PD1dXV7Zv386MGTNYt24d/v7+jX7t/PxSh7q8LTjYh9zcEluX0aTUk/1wxL7sraeKGjM5JVXkllaTU1pV7/Nf/i2rMvO3xwbQ3tP1kl/fZHIiMNC7weU2GScxYcIEXnzxRQoLC+sFQHBwcN3ngwcPJjw8nJSUFAYMGGCLMkVEmlVFjZnMs5WcKa4ku6SKnNJqcn8OgezSKnJLqyitMp/3PG93Z4K93QnxdqNzoD/hvh50C/WhqqSyyWu0SkiUlZVRXFxMeHg4AJs3b8bX1xc/P79662VnZxMaGgrAoUOHyMjIoFOnTtYoUUSkydWaLWSVVJF5tvLcR/G5fzN+/rqgvKbe+iYnCPRyI9jbnQ7+begf5UeItxshPu4Ee7v9HAzueLo5n7etth6u5NprSFRUVDB9+nQqKiowmUz4+vqydOlSnJycSEhIYNq0acTExDBv3jwOHDiAyWTC1dWVN954o97ehYhIS1NcWcOpggoyzlb8Kwx+/sguqcL8q6Pfzk4Q2taDCF8PbuwcSISvR91HqI87gV5uuJha1kU9Djd9qc5JtHzqyX44Yl+X05PFMMgpqeJkQTmnCio4XVDOyfxyThWUn7c3EOjlRkRbDyJ83Wn3cwC0821DhK8HIT7uzRYCl/u9apHnJEREWqLqWgtpRRWcKjgXACfzyzldcO7rylpL3XptPVzoGODJDZ0D6BjgSYcAT6L82hDe1h0P1/MPBdkzhYSItDoWwyC9qJKU3FKO5pZxLLeMUwXlZBRV1Ds8FN7WnQ4BntwWGU7HgDZ0DPSkY4An/m1cW81YL4WEiDi0supajuWWcTS3jJTcUk4WVnL4THHdnoGzE7QP8CQ62IsR3YPpGOBJpwBP2ge0oY2D7RVcDoWEiDgEwzDILK4kJaeMlNwyjuaWkpJbRsbZf13x4+PuQs+ItsTHhBEd7E10iBedAr1wd9G9ThuikBARu/PL4aKDWSUcyi7hUFYJR3PLKKs+N6bACYjyb0OPUG/G9wqja7AX0cFehPq4ExLS1uFOxjcnhYSItGiGYZBdUsXBrBIOZJWeC4XskrpBZu4uJqKDvbmlRwjdQryJDvaiS5CXDhU1EYWEiLQo+WXVHMwq+XkvoZSDWSUUVpy7zNTF5ES3YC9Gdg+hZ5g3PUJ96Bzk1eLGFjgShYSI2Ex1rYWDWSX8mHGWAz8HQ05pNXBu9HGnQE8Gdw6gZ5gPPcN86Bqk8wfWppAQEauprDGz/0wJP6QXsSf9LD+dKaHq56uM2vu3oV+k77lACPUhOsT7grefEOtSSIhIsymvNrMv8yw/pJ/lh7Rzewu1FgOTE0QHe/O7PuH0a+dL30hf/Npc+h1MpfkpJESkyZRU1vJjxln2pJ8LhsPZJZiNc2MReoT5cG9sO/pF+tK3nS/e7nr7sQf6LonIZauutfBDehHfnipkd9pZjuaUYgCuzk5cHebDgwOi6BfpS+8IXx06slMKCRG5JDklVWw/WcD2EwXsTC2kosaCu4uJmHAfEgZ1oF+kL73CfRzuHkatlUJCRC7KbDE4mFXCtp+D4UhOKQBhPu6M6RnKDZ0D6B/lp1BwUAoJETlPSWUt350uZNfm42w+nENRRQ0mJ+gd0ZbJN3Tkhi6BdAn0bDU3uWvNFBIigmEYnCqoYNuJfLadKGBvxlnMBvh5unJdR39u6BTAdR398dUVSK2OQkKklTIMgyM5pWw6ksuXR/PqboTXLdiLSddGcUPnAIb2bkdBfqmNKxVbUkiItDLH88rYdCSXTUdySS2swNnkxMAOfjxwbSTXdwogrK1H3brOut1Fq6eQEGkFUgsr2HQkh01HcjmeV47JCWKj/Li/fyRDuwVpIJs0SCEh4qDOFFeSfCSXjYdzOfzzFUl927XlT3FdGRYdRKCXm40rFHugkBBxIHmlVSQfzWPj4Vx+OlMMwNVhPvx+SGeGRQfVO5Qk0hgKCRE7V1xZc26P4UguP6SdxeDcyefJN3RkePdgIv3a2LpEsWMKCRE7ZLYYfJ9ayJr92fzzWB7VZoOOAW1IGNTh3DzNgZ62LlEchEJCxI6kF1Ww5kA2a/dnkVNaja+HC7f1Dmfc1WFEh3hpcJs0OYWESAtXUWPmy6O5rN6fzZ70s5ic4LqO/jx9cxdu6hKImybhkWakkBBpgQzDYF9mMWv2Z7PpSC7lNWai/Dx46oaO3NozlBAfd1uXKK2EQkKkBcktrSLpQDZrDmSTWlhBG1cTI7oHM+7qMPq0a6vDSWJ1CgkRG6s1W/jqeD5r9mfz7akCLAb0i/TloQFRDIsO1jwMYlMKCREbKa6sYdW+LD7dk0FOaTUh3m48NCCKsVeHEeWvy1alZVBIiFhZWmEF//NDBmsOZFFRY6F/ez+eHd6N6zsF6F5J0uIoJESswDAM9mSc5ZNdGXx9PB9nkxOjeoRw7zXtiA7xtnV5Ig1SSIg0o1qzheSjeXyyO51D2aX4erjw8HXtmdgnnCBvXaEkLZ9CQqQZFFfWsHJfFv/78/mGDv5tmDW8K2N6hmqaT7ErVguJp556ivT0dEwmE56enrzwwgv06NGj3jpms5lXX32VrVu34uTkxOOPP87EiROtVaLIFUv95XzD/iwqay1c296P50ZEM6iTPyZdvip2yGohkZiYiI+PDwDJyck899xzrFy5st46a9asITU1lY0bN1JUVMSECRMYNGgQkZGR1ipT5JIZhsF3J/JZ8mUKW4/n4+LsxKirQrhH5xvEAVgtJH4JCIDS0tILDgpat24dEydOxGQyERAQwPDhw1m/fj2PPfaYtcoUuSR7M86yeNsp9qSfxdfDhUeua88dfSMI0lwN4iCsek7i+eefZ/v27RiGwbvvvnve8jNnzhAREVH3dXh4OFlZWdYsUaRRjuaU8t/bT7HtRAEBnq78eVxPhnXy1/kGcThWDYnXXnsNgFWrVvHGG2+wbNmyJt9GYKDj7d4HB/v89kp2xl57OplXxrxNR1mzN5O2Hi7MHN2dh67viKeb414DYq/fq4txxJ6gefqyyU/2hAkTePHFFyksLMTf37/u8fDwcDIzM+nduzdw/p5FY+Tnl2KxGE1ary0FB/uQm1ti6zKalD32lF1SxbvfnmbN/ixcnU08PDCK+/tH0tbDlbKzFXjaYU+NYY/fq9/iiD3B5fdlMjld9I9rq4REWVkZxcXFhIeHA7B582Z8fX3x8/Ort97o0aNZsWIFI0eOpKioiOTkZP7+979bo0SRCyosr+b/7Uzj/37MxADu6BvBwwPba35oaTWsEhIVFRVMnz6diooKTCYTvr6+LF26FCcnJxISEpg2bRoxMTHEx8ezd+9eRo4cCcDkyZOJioqyRoki9ZRW1fL3Xel8sjuDylozt/YMJeH6DoRrjmhpZZwMw3CcYzPocJM9aMk9VdaYWfFjJh/uTONsZS3DooN44vqOdPqN6UBbck9XwhH7csSewM4PN4m0dLVmC6v3Z/Hud6nkllZzXUd/nrqhIz1CHfMEp0hjKSSk1fv6eD7z/3mc9KJKeke05ZUxVxEb5ffbTxRpBRQS0mplnq3kr1uO8/XxfDoHejL/tqsZ3ClAs7+J/IpCQlqdGrOFj3el8953qZicYNpNnbjnmna4OJtsXZpIi6OQkFZlV2oRiV+mcKqggqHdgvjDzZ0J0xVLIg1SSEirkFdWzZtfnWD9oRza+Xqw4LZeDO4cYOuyRFo8hYQ4NLPF4LO9mSzZdopqs4XHrmvPgwOidI8lkUZSSIjDOnCmmLnJxzicU8qA9n7MHNaVDgEXH+8gIvUpJMThFFfWsGTbKT7fe4ZALzfmjO3B8OggXbUkchkUEuIwDMNg3cEc3vzqBGcra7j7mnY8fn0HvN31Yy5yufTbIw7heF4Zickp7MkoJia8LQuHx9Bds8KJXDGFhNi1GrOF979L5YOdaXi7OfP8iG6MjwnTfNIiTUQhIXbrWG4ZL31xmKO5ZdzSI4Q/3NwFP09XW5cl4lAUEmJ3ai0Gf/s+jXe+OU1bDxf+Et+TIV2DbF2WiENSSIhdOZVfzp/XH+FAVgnDo4N4Zlg37T2INCOFhNgFi2HwPz9ksGTbKTxcTLx261WMvCrE1mWJODyFhLR46UUVzF5/hD0ZxdzUJZBZI7oRpOlDRaxCISEtlsUw+GzvGd766gQuzk68NDqaW3uGalCciBUpJKRFyiqu5JUNR9mZWsR1Hfz5z1HRhPq427oskVZHISEtimEYrNmfzbx/HscwYNaIbtwWE6a9BxEbUUhIi5FXWsVrm1LYdqKA2ChfXhgVTTvfNrYuS6RVU0iIzRmGwYbDufzX5mNU1Vr449Au3NkvQqOmRVoAhYTYVFWthbnJKaw9kE1MeFteGh2t23mLtCAKCbGZrOJKZq4+yKHsUh67rj2PDeqAs0l7DyItiUJCbGJXahGz1h6ixmzhL/FXM6RroK1LEpELUEiIVRmGwbtbT/D6ukO09/fkjfiedNThJZEWSyEhVlNZY+bVjUfZcDiXod2CeGl0NF5u+hEUacn0GypWkV5UwczVBzmWW8afRnVn4tUhGvsgYgcUEtLsvjtVwPNJhzEMWHB7L+IHdCA3t8TWZYlIIygkpNkYhsGHO9NYsu0UXYK8+K/4nkT6aXCciD1RSEizKKuu5ZUNR/nyaB4juwfzn6OiaePqbOuyROQSKSSkyZ0uKOdPqw9yuqCc3w/pzL2x7XT+QcROKSSkSW09ns8L6w7j6mxi0R0xXNve39YlicgVUEhIk7AYBu99m8o7357mqhBv3ojvSXhbD1uXJSJXyCohUVhYyMyZM0lNTcXNzY0OHTowe/ZsAgIC6q23cOFCPvnkE0JCzk1Lec011/DSSy9Zo0S5AqVVtby47jBbTxRwa88Qnh3eDQ+dfxBxCFYJCScnJx577DEGDhwIQGJiIn/5y1+YM2fOeetOmDCBZ555xhplSRMoqqhh2mc/cTS3jD/FdWFi3widfxBxICZrbMTPz68uIAD69u1LZmamNTYtzSivrJon/3cvx/PK+K/xPbmzn05Qizgaq5+TsFgsLF++nLi4uAsuT0pKYtu2bQQHBzN16lT69et3Sa8fGOjdFGW2KMHBPrYu4TyZRRU89X+7yTpbxQcPD2Bw16BLen5L7OlKOWJP4Jh9OWJP0Dx9ORmGYTT5q17Eyy+/THZ2NosWLcJkqr8jk5ubi5+fH66urmzfvp0ZM2awbt06/P0bf4VMfn4pFotVW2pWwcE+LW50clphBZP/bx/FlbW8eXsv+rTzvaTnt8SerpQj9gSO2Zcj9gSX35fJ5HTRP66tcrjpF4mJiZw+fZoFCxacFxAAwcHBuLq6AjB48GDCw8NJSUmxZonyG07kl/H4p3sprzaz9M7elxwQImJfrBYS8+fPZ//+/SxevBg3N7cLrpOdnV33+aFDh8jIyKBTp07WKlF+w+HsEp74dB8G8PZdfbgq1DF32UXkXxp9TsIwDFasWMHatWspLCxkzZo1fP/99+Tm5jJmzJiLPjclJYWlS5fSsWNH7r77bgAiIyNZvHgxCQkJTJs2jZiYGObNm8eBAwcwmUy4urryxhtvEBwcfGUdSpPYl1nM9M9/wtvNhSUTexPlr3swibQGjQ6JN998k2+++YYHH3ywbuxCWFgYr7/++m+GRLdu3Thy5MgFly1btqzu88TExMaWI1b0fWohf1x1gCAvN5ZM7E2YBsmJtBqNPty0cuVKli5dyq233lp3mWNkZCRpaWnNVpzY3rYT+fz+8/2Et/Xgnbv6KCBEWplG70mYzWa8vLwA6kKirKwMT09NPemoko/k8p/rDhMd7MVbv4vBr42rrUsSEStr9J7EkCFDeP3116murgbOnaN48803GTp0aLMVJ7az9kAWzycdoleYD0sm9lZAiLRSjQ6JWbNmkZOTQ2xsLCUlJfTr14/MzExmzJjRnPWJDaz4MZOX1x+lf5QfC++Iwdtd94EUaa0a/dvv7e3NkiVLyM/PJyMjg/DwcF155ID+9n0ab319kpu6BDJnbA/cXaw6lEZEWphGh0RBQQHu7u4EBgbi5+fHqlWrcHZ2Zvz48RccGCf2xTAM3vnmNO9+l8qI7sHMvqU7Ls76voq0do1+F3jiiSc4ffo0cG5g3Pvvv88HH3zA3Llzm604sQ7DMFjw1Qne/S6V8b1CeWXMVQoIEQEuISROnTpFjx49AFi9ejXLli3jww8/ZN26dc1WnFjH4m2n+GR3Bnf1i+D5kdE4m3QnVxE5p9GHm0wmEzU1NZw8eRIfHx8iIiKwWCyUlZU1Z33SzD7fm8mHO9P4XZ9w/ji0i271LSL1NDokbrrpJqZPn05RUVHdCOtjx44RGhrabMVJ89p+ooDEL49xQ+cAZsR1VUCIyHkaHRKvvfYaK1euxMXFhfj4eODctKRTp05ttuKk+RzKLmHW2oNEB3vz2q09cNEhJhG5gEaHhJubG3fddVe9x34925zYjzPFlTy98gB+bVyZf3svPN00H7WIXFijQ6KkpISPPvqIQ4cOUV5eXm/Z+++/3+SFSfMoqaxl+uf7qao1s2RiX4K8LnzbdhERuISQmD59OmazmREjRuDu7t6cNUkzqa618KfVB0grrGDRHTF0DvSydUki0sI1OiR+/PFHduzYUTdznNgXwzB4deNRdqedZfaY7sRG+dm6JBGxA40eJxEbG8vx48ebsxZpRku/Oc0Xh3J46oaO3NJDV6SJSOM0ek9i7ty5JCQk0KdPHwIDA+stmzJlSpMXJk1n1b4zvP9dKvExYTw0IMrW5YiIHWl0SMyfP5+srCwiIyMpLS2te1zX1rds354qYG5yCoM6+vPsMI2FEJFL0+iQSEpKYsOGDYSEhDRnPdKEjuSU8uzqQ3QJ8uL1cT10PyYRuWSNfteIiorCxUXzCtiLrOJKnl65Hx8PFxbc3gsvN33vROTSNfqdIz4+nqeeeor777//vHMSgwYNavLC5PKVVtXy+5X7Ka828+49fQn21iXLInJ5Gh0SH330Ec7OzsybN6/e405OTnz55ZdNXphcnhqzhWdWH+RUQQVv3d6LrkEaCyEil69RIWE2mykoKGD37t24uWmEbktlGAavbUphZ2oRL42OZkAHf1uXJCJ2rlHnJJydnenUqROFhYXNXY9cgXe/TSXpQDaPX9+BsVeH2bocEXEAjT7cNG7cOJ588kkeeOABwsLqvwHpnITtrdmfxTvfnmbs1aE8dl17W5cjIg6i0SGxfPlyABYuXFjvcZ2TsL1dqUW8timFgR38eH5EN42FEJEm0+iQ2Lx5c3PWIZcpq7iS59Yeor1fG+aO66mxECLSpPSOYseqay3MWnuIarOFN+J74u2usRAi0rQUEnZs3j+Ps/9MCS+N7k7HAE9blyMiDkghYafW7M/is71neODaKIZ2C7J1OSLioBQSduhIdimJXx6jf3s//uOGjrYuR0QcmELCzpytqGHm6gP4ergw59arcDHpSiYRaT4602lHLIbBC+sOk1tWzbK7+uDvqdHvItK8rBIShYWFzJw5k9TUVNzc3OjQoQOzZ88mICCg3npms5lXX32VrVu34uTkxOOPP87EiROtUaJdWPbNab49Vcis4V25OrytrcsRkVbAKoebnJyceOyxx9iwYQNr1qwhKiqKv/zlL+ett2bNGlJTU9m4cSOffvopCxcuJD093Roltnhbj+fz7nepjLs6lNt6h9u6HBFpJawSEn5+fgwcOLDu6759+5KZmXneeuvWrWPixImYTCYCAgIYPnw469evt0aJLVp6UQUvfXGE7iHezNTsciJiRVY/cW2xWFi+fDlxcXHnLTtz5gwRERF1X4eHh5OVlWXN8lqcimozM1cfxMkJEsf3wMPV2dYliUgrYvUT16+88gqenp7cf//9zfL6gYHezfK6tmAYBn/4370cyyvjg4eupW9Xx5k6NjjYx9YlNDlH7Akcsy9H7Amapy+rhkRiYiKnT59m6dKlmEzn78SEh4eTmZlJ7969gfP3LBojP78Ui8VoknptbcWPmazck8ET13fg6oA25OaW2LqkJhEc7OMwvfzCEXsCx+zLEXuCy+/LZHK66B/XVjvcNH/+fPbv38/ixYsbnLho9OjRrFixAovFQkFBAcnJyYwaNcpaJbYo+zKLmbflOMOuCuER3fpbRGzEKiGRkpLC0qVLycnJ4e677yY+Pp7JkycDkJCQwE8//QScm0c7MjKSkSNHcueddzJ58mSioqKsUWKLkl9WzbNrDhLW1p15d/XFpBPVImIjToZhOMaxmZ/Z++GmWovB5BX7OJBVwvv39GXw1eEOt2vsiLv7jtgTOGZfjtgTOMDhJmmcRV+f5If0szw3ohvRIY5zEl5E7JNCogVJPpLL33enc2ffCMb0DLV1OSIiComW4mR+ObM3HCEmvC2/v7mzrcsREQEUEi1CaVUtf/rHAdq4OjN3XA9cNQWpiLQQugtsC/DW1ydIK6pgycTehPi427ocEZE6+pPVxvZmnGXlvizuuSaS2Cg/W5cjIlKPQsKGaswWXtuUQpiPO08M7mDrckREzqOQsKGPd6VzMr+cmcO60kY37hORFkghYSNphRW8910qw6KDuLFLoK3LERG5IIWEDRiGQeKXKbiYnPjj0C62LkdEpEEKCRtYfziHHaeLmHxjJ4K9dTWTiLRcCgkrO1tRw/wtJ+gV7sPtmoZURFo4jZOwsoVbT1JcWcOi4TE4m3R3VxFp2bQnYUV70s/yj5+yuDc2UjfvExG7oJCwkupaC69vSiGirTsJ12tMhIjYB4WElfxtVxonC8YltBoAAA3fSURBVMqZObybxkSIiN1QSFhBamEF73+XyvDoYAZ3CrB1OSIijaaQaGaGYTA3OQU3FxN/HKpbgIuIfVFINLMvDuXwfWoRU27sRJDGRIiInVFINKOiihrm//MEMeE+3KYxESJihxQSzWjh1ycoqarluRHRmJw0JkJE7I9CopnsTiti9f5s7ouNpGuwl63LERG5LAqJZlA3JsLXg4RB7W1djojIZVNINIMPv0/jdGEFzwzriofGRIiIHVNINLHTBeV8sCOVkd2DuV5jIkTEzikkmtAvYyLcXUw8rXkiRMQBKCSaUNLBbHalnWXqjZ0I8nKzdTkiIldMIdFEisprWPDPE/SOaMsEjYkQEQehkGgib359gtJqM7NGdNOYCBFxGAqJJvBDehFrD2QzqX8kXYM0JkJEHIdCogl8sCONIC83Hr1OYyJExLEoJK7QifwyvjtVyMS+ERoTISIORyFxhZbvzsDdxcTtOlktIg7IKiGRmJhIXFwc3bt35+jRoxdcZ+HChQwaNIj4+Hji4+N5+eWXrVHaFSkqr+GLQznc0iMEP09XW5cjItLkXKyxkWHDhvHAAw9w3333XXS9CRMm8Mwzz1ijpCbx+b4zVNVauCe2na1LERFpFlYJif79+1tjM1ZVY7aw4sdMruvoT+dAXdEkIo6pRZ2TSEpKYty4cTzyyCPs2bPH1uVc1KYjueSVVXPPNdqLEBHH5WQYhmGtjcXFxbF06VKio6PPW5abm4ufnx+urq5s376dGTNmsG7dOvz9/a1VXqMZhsHYhduoqrWw6embcNLgORFxUFY53NQYwcHBdZ8PHjyY8PBwUlJSGDBgwCW9Tn5+KRZL8+beD+lFHMgsZtaIbuTllTbrtoKDfcjNLWnWbViberIfjtiXI/YEl9+XyeREYKB3w8uvpKimlJ2dXff5oUOHyMjIoFOnTjasqGHLd2fg6+HCmB4hti5FRKRZWWVP4tVXX2Xjxo3k5eXx8MMP4+fnR1JSEgkJCUybNo2YmBjmzZvHgQMHMJlMuLq68sYbb9Tbu2gp0osq+OpYPg8NjNLgORFxeFY9J2ENzX246a9bjrPix0zWJAwg2Nu92bbzC0fcNVZP9sMR+3LEnqAVHG6yB6VVtaz+KYsR3YOtEhAiIramkLgEq/dnUV5j5l4NnhORVkIh0Uhmi8GnP2TQr11beoT62LocERGrUEg00lfH8sgsruLu2EhblyIiYjUKiUb6ZHcGEb4eDOkSaOtSRESsRiHRCAeyStibWcxd/SJwNml0tYi0HgqJRli+Ox0vN2fG9wqzdSkiIlalkPgNOSVVJB/NY3yvMLzdW8xdTERErEIh8RtW/JiJYRjcdU2ErUsREbE6hcRFVNaYWbnvDDd1CaSdbxtblyMiYnUKiYtIOpjN2cpa7tVlryLSSikkGmAxDP7nhwx6hHrTt11bW5cjImITCokGfHuqkFMFFdx9TTtNKiQirZZCogHLd6cT5OXGiO4t73blIiLWopC4gON5Zew4XcSd/SJwddZ/kYi0XnoHvIDlP2Tg7mLitphwW5ciImJTCol/U1hezRcHsxnTMwQ/T1dblyMiYlMKiX/z2d4zVJsN7r5Gc0aIiCgkfqW61sL/7T3DdR396RzoZetyRERsTiHxK5uO5JJfVq2Z50REfqaQ+JlhGHyyO51OAZ5c18Hf1uWIiLQIComf/ZB+lqO5Zdwdq8FzIiK/UEj8bPnuDHw9XBjTI8TWpYiItBgKCSCtsIKvj+fzuz7heLg627ocEZEWQyEBHM0txdPNmTv6as4IEZFf01RrwLDoYK7r6I+Xm/47RER+TXsSP1NAiIicTyEhIiINUkiIiEiDFBIiItIghYSIiDRIISEiIg1SSIiISIMc7rpPk8nx7ruknuyDI/YEjtmXI/YEl9fXbz3HyTAM43ILEhERx6bDTSIi0iCFhIiINEghISIiDVJIiIhIgxQSIiLSIIWEiIg0SCEhIiINUkiIiEiDFBIiItIghUQLUFhYSEJCAqNGjWLcuHFMmTKFgoICAE6ePMldd93FqFGjuOuuuzh16pRti71MixYtonv37hw9ehSw776qqqp46aWXGDlyJOPGjeOFF14A7LsngC1btjBhwgTi4+MZN24cGzduBOyrr8TEROLi4ur9rMHFe2jp/V2op4u9Z0AT92SIzRUWFhrfffdd3ddz5841Zs2aZRiGYUyaNMlYtWqVYRiGsWrVKmPSpEk2qfFK7N+/33j00UeNm2++2Thy5IhhGPbd1yuvvGK89tprhsViMQzDMHJzcw3DsO+eLBaL0b9//7rvz6FDh4y+ffsaZrPZrvr6/vvvjczMTGPo0KF1vRjGxb83Lb2/C/V0sfcMw2janhQSLdD69euNBx980MjLyzNiY2ON2tpawzAMo7a21oiNjTXy8/NtXGHjVVVVGXfeeaeRmppa90Nuz32VlpYasbGxRmlpab3H7bknwzgXEgMGDDB27dplGIZh7Ny50xg5cqTd9vXrN9SL9WBP/f178P3aL+8ZhtH0P4sOdxdYe2exWFi+fDlxcXGcOXOG0NBQnJ2dAXB2diYkJIQzZ84QEBBg40ob580332T8+PFERUXVPWbPfaWlpeHn58eiRYvYsWMHXl5eTJ8+HQ8PD7vtCcDJyYkFCxbw1FNP4enpSVlZGW+//bZdf69+cbEeDMOw+/5+/Z4BTf/7pXMSLcwrr7yCp6cn999/v61LuWJ79uzhp59+4t5777V1KU2mtraWtLQ0evbsyeeff86MGTOYOnUq5eXlti7titTW1vL222+zZMkStmzZwn//93/z9NNP231frUFzv2coJFqQxMRETp8+zYIFCzCZTISHh5OdnY3ZbAbAbDaTk5NDeHi4jSttnO+//54TJ04wbNgw4uLiyMrK4tFHHyU1NdVu+4qIiMDFxYWxY8cC0KdPH/z9/fHw8LDbngAOHTpETk4OsbGxAMTGxtKmTRvc3d3tui/gor9H9v479u/vGXDxfi+HQqKFmD9/Pvv372fx4sW4ubkBEBgYSI8ePVi7di0Aa9eupUePHnazG/z444+zbds2Nm/ezObNmwkLC+O9995jzJgxdttXQEAAAwcOZPv27cC5q0jy8/Pp2LGj3fYEEBYWRlZWFidOnADg+PHj5OXl0aFDB7vuCy7+e2TPv2MXes+Apn/f0KRDLUBKSgpjx46lY8eOeHh4ABAZGcnixYs5fvw4zz77LMXFxbRt25bExEQ6d+5s44ovT1xcHEuXLiU6Otqu+0pLS+O5556jqKgIFxcXfv/73zNkyBC77glg9erVLFu2DCenczOVTZs2jeHDh9tVX6+++iobN24kLy8Pf39//Pz8SEpKumgPLb2/C/W0YMGCBt8zoGl7UkiIiEiDdLhJREQapJAQEZEGKSRERKRBCgkREWmQQkJERBqkkBARkQYpJKRViouL45tvvmnS1/z888+55557Lvl5u3btYtSoUU1ai0hTUUiIWFn37t05ffp03df9+/dnw4YNNqxIpGEKCRERaZBCQlo1i8XCO++8w/Dhwxk4cCDTp0+nqKiobvm0adMYPHgwsbGx3HfffaSkpNQtKyws5Mknn+Saa67hjjvuIDU19Te3d9999wEQHx9Pv379WLduHTt27OCmm26qWycuLo53332XcePG0bdvX5577jny8vJ47LHH6NevHw899BBnz56tW//HH3/k7rvvpn///owfP54dO3Y0xX+NCKCQkFbuo48+Ijk5mY8//pitW7fi6+vL7Nmz65bfdNNNbNiwgW+//ZaePXsyY8aMumWzZ8/G3d2dbdu2MWfOHD777LPf3N7f//53AP7xj3+wZ88exowZc8H1Nm7cyAcffMCGDRvYsmULCQkJ/OEPf2DHjh1YLBb+9re/AZCdnc0TTzzBf/zHf7Bz506eeeYZpk2bVm8qS5EroZCQVu3TTz/l6aefJiwsDDc3N6ZMmcKGDRuora0F4I477sDb2xs3NzemTp3K4cOHKSkpwWw2s3HjRqZNm4anpyfR0dHcdtttTVbX/fffT1BQEKGhofTv35/evXvTs2dP3NzcGDFiBAcPHgTOhc1NN93EkCFDMJlMDB48mF69evHVV181WS3SumlmOmnVMjMzmTx5ct29+AFMJhP5+fkEBQUxf/581q9fT0FBQd06hYWFVFZWUltbW+8e/REREU1WV1BQUN3n7u7u9b728PComwwoMzOT9evXs2XLlrrltbW1DBw4sMlqkdZNISGtWlhYGHPmzKmbbOfXVq1axZdffskHH3xAZGQkJSUlXHvttRiGQUBAAC4uLpw5c4YuXboA56aNtLbw8HDi4+N59dVXrb5taR10uElatXvuuYcFCxaQkZEBQEFBAcnJyQCUlZXh5uaGv78/FRUVzJs3r+55zs7OjBgxgkWLFlFRUcGxY8dYuXJlo7YZFBREWlpak9Q/fvx4tmzZwtatWzGbzVRVVbFjxw6ysrKa5PVFFBLSqj3wwAPExcXxyCOP0K9fP+6880727dsHwIQJE4iIiODGG2/k1ltvpW/fvvWe++KLL1JeXs7gwYN59tlnuf322xu1zSlTpvDss8/Sv39/1q1bd0X1h4eHs2TJEt5++20GDRrEkCFDeO+997BYLFf0uiK/0KRDIiLSIO1JiIhIg3TiWqSJ7dq1i4SEhAsu27Nnj5WrEbkyOtwkIiIN0uEmERFpkEJCREQapJAQEZEGKSRERKRBCgkREWnQ/wev3vKZGHNHIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rmses_t850_6h_iter.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'rmse' (lead_time: 2)>\n",
       "dask.array<getitem, shape=(2,), dtype=float64, chunksize=(1,)>\n",
       "Coordinates:\n",
       "  * lead_time  (lead_time) int64 72 120"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rmses_z500_6h_iter.sel(lead_time=[3*24, 5*24])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'rmse' (lead_time: 2)>\n",
       "dask.array<getitem, shape=(2,), dtype=float64, chunksize=(1,)>\n",
       "Coordinates:\n",
       "    level      int32 850\n",
       "  * lead_time  (lead_time) int64 72 120"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rmses_t850_6h_iter.sel(lead_time=[3*24, 5*24])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# The End"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
